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One  of  the  most  important  problems  in  astronomy  today  is  the  development  of 
a comprehensive  theory  for  the  origins  of  stars  and  planetary  systems.  In  order  to  un- 
derstand star  formation,  we  must  understand  the  interstellar  medium  (ISM)  in  which 
the  process  occurs.  The  evolution  of  the  ISM  is  driven  by  a variety  of  phenomena, 
including  gravity,  turbulence,  shearing  flows,  magnetic  fields,  and  the  thermal  prop- 
erties of  the  gas.  As  a result,  the  ISM  is  susceptible  to  a wide  range  of  hydrodynamic 
instabilities.  In  this  dissertation,  two  very  different  types  of  instabilities  are  studied 
using  computational  hydrodynamics  simulations.  The  first  is  a recently  discovered 
instability  which  acts  upon  an  interface  of  discontinuous  density.  Theory  predicts 
that  this  self-gravity  driven  interfacial  instability  persists  in  the  static  limit  and  in 
the  absence  of  a constant  background  acceleration.  Disturbances  to  a density  inter- 
face are  found  to  grow  on  a time-scale  of  the  order  of  the  free-fall  time,  even  when  the 
perturbation  wavelength  is  much  less  than  the  Jeans  length.  Here,  I present  the  first 


numerical  simulations  of  this  instability.  The  theoretical  growth  rate  is  confirmed  and 
the  nonlinear  morphology  displayed.  The  self-gravity  interfacial  instability  is  shown 
to  be  fundamentally  different  from  the  Rayleigh- Taylor  instability,  although  both  ex- 
hibit similar  morphologies  under  the  condition  of  a high  density  contrast  such  as  is 
commonly  found  in  the  ISM.  Such  instabilities  are  a possible  mechanism  by  which 
observed  features,  such  as  the  pillars  of  gas  seen  near  the  boundaries  of  interstellar 
clouds,  are  formed. 

The  second  instability  considered  in  this  study  is  a dynamic  instability  which 
occurs  when  supersonic  gas  streams  collide.  It  is  driven  by  a combination  of  ram 
pressure  and  strong  cooling.  I examine  the  stability  of  a cold,  dense  slab  formed  by 
the  head-on  collision  of  supersonic  gas  streams.  Special  emphasis  is  placed  on  the 
effects  of  altering  the  cooling  rate.  A series  of  models  are  performed  in  cylindrical 
geometry  to  determine  whether  the  instability  of  the  slab  acts  to  increase  or  decrease 
the  likelihood  of  forming  collapsing  regions  capable  of  undergoing  star  formation. 
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CHAPTER  1 
INTRODUCTION 


One  of  the  most  important  problems  in  astronomy  today  is  the  development 
of  a comprehensive  theory  for  the  origins  of  stars  and  planetary  systems.  This  is  a 
problem  of  great  complexity  that  involves  changes  over  many  orders  of  magnitude  in 
size,  density,  and  temperature.  Advances  in  observational  and  computational  tools 
are  fueling  progress,  but  many  questions  remain  unanswered.  In  order  to  understand 
star  formation,  we  must  understand  the  environment  in  which  the  process  occurs- 
molecular  clouds.  The  morphology  and  evolution  of  molecular  clouds  are  driven 
by  a variety  of  physical  phenomena,  some  of  which  aid  the  star  forming  process 
while  others  oppose  it.  This  dissertation  describes  a theoretical  study  that  uses 
computational  techniques  to  investigate  hydrodynamic  instabilities  which  may  trigger 
the  formation  of  self-gravitating  regions  conducive  to  star  formation.  In  this  chapter,  I 
summarize  current  theories  of  molecular  cloud  morphology  and  star  formation.  I also 
give  a brief  history  detailing  the  use  of  numerical  simulations  to  study  hydrodynamic 
flows  in  the  interstellar  medium  (ISM).  The  hydrodynamic  instabilities  that  are  the 
focus  of  this  dissertation  are  introduced.  This  chapter  concludes  with  an  outline  for 
the  remainder  of  the  dissertation. 

1.1  Molecular  Clouds  and  Star  Formation 

Because  all  present  day  star  formation  is  observed  to  take  place  inside  giant 
molecular  clouds  (GMCs),  knowing  how  molecular  clouds  form  and  evolve  is  crucial 
for  understanding  the  initial  conditions  for  star  formation.  Observational  work  using 
radio  and  infrared  telescopes  has  produced  a wealth  of  information  about  GMCs.  A 
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typical  GMC  is  gravitationally  bound,  magnetized,  and  compressible.  GMCs  in  the 
solar  neighborhood  are  observed  to  have  masses  of  10^“®  Mq  and  a mean  diameter 
of  45  pc  (Blitz  1993).  Molecular  clouds  are  observed  to  be  clumpy.  The  clumps  vary 
in  geometry  from  spherical  to  highly  filamentary.  Molecular  linewidths  are  wider 
than  thermal  linewidths  inferred  from  excitation  temperatures,  indicating  a turbulent 
medium.  The  observed  structures  of  molecular  clouds  are  self-similar  on  different 
length  scales  and  obey  power-law  scaling  relations  (Larson  1981).  The  relations 
among  the  internal  velocity  dispersion  Au,  mean  gas  particle  density  n,  cloud  radius 
R,  and  cloud  mass  M are  (from  Falgarone  and  Puget  1986) 

( (1,1) 

\0.36kms~^J  \570  cm~^/  \pc  J \100  MqJ 

Recent  evidence  reveals  that  the  self-similarity  does  not  extend  to  small  (<  0.2pc) 
scales  (Williams  1999),  indicating  that  once  a core  begins  to  collapse  due  to  self- 
gravity it  takes  on  characteristics  very  different  from  the  parent  cloud. 

Over  the  last  century,  star  formation  has  been  the  subject  of  extensive  scien- 
tific study.  According  to  the  classic  theory  of  Jeans  (1929),  a static,  homogeneous, 
spherical  cloud  will  collapse  if  its  meiss  exceeds  a critical  value,  Mj,  given  by 

= r'’  (^)  ’ ■ 

In  the  above,  p is  the  gas  volume  density,  c the  sound  speed,  and  G the  gravitational 
constant.  The  time  scale  for  collapse  is  the  free-fall  time,  iff  = {Zir /2>2Gp)^^^^\ 
For  densities  typical  of  molecular  clouds,  a free-fall  time  of  tjf  ~10®?/r  is  expected. 
Although  no  current  star  formation  is  observed  outside  of  molecular  clouds,  there 
are  clouds  in  which  no  star  formation  is  observed.  This  presents  two  possibilities: 
either  there  are  some  clouds  in  which  conditions  are  intrinsically  unsuitable  for  star 
formation,  or  there  are  clouds  in  which  star  formation  has  not  yet  begun.  This  second 
possibility  implies  that  molecular  clouds  have  lifetimes  that  exceed  the  time  needed 
for  star  formation.  If  clouds  are  to  have  lifetimes  longer  than  their  free-fall  times. 
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some  means  of  support  against  gravitational  collapse  is  needed.  One  possible  support 


clouds,  determined  by  measuring  the  Zeeman  effect,  show  a range  from  10  - 3100  jiG 
(e.g.  Crutcher  1999).  Magnetic  fields  can  support  a cloud  provided  the  cloud  mass 
does  not  exceed  a critical  value,  M^it,  such  that  (Lada  and  Shu  1990) 


According  to  the  popular  bimodal  theory  of  star  formation,  this  critical  mass  di- 
vides two  regimes  of  star  formation.  Clouds  that  exceed  the  critical  mass  are  called 
supercritical,  and  clouds  with  a mass  below  the  critical  mass  are  called  subcritical. 
In  a supercritical  cloud,  magnetic  support  can  not  prevent  collapse.  The  cloud  will 
contract  and  fragment  into  smaller  clouds  that  may  or  may  not  be  supported  by 
magnetic  fields  themselves.  For  cases  in  which  the  magnetic  field  strength  is  low,  a 
large  cloud  that  contracts  and  fragments  may  ultimately  form  a cluster  of  stars.  If 
little  or  no  fragmentation  occurs,  the  process  produces  high  mass  stars.  In  subcritical 
clouds,  magnetic  fields  slow  the  collapse.  Ions  are  coupled  to  the  lines,  but  neutral 
molecular  gas  is  free  to  diffuse  across  magnetic  field  lines.  When  neutral  molecules 
attempt  to  cross  a field  line  due  to  gravitational  contraction,  their  motion  is  retarded 
by  collisional  interaction  with  the  ions.  This  process  of  ambipolar  diffusion  slows 
the  cloud  contraction  so  that  the  time  needed  for  collapse  is  increased  by  an  order 
of  magnitude.  The  typical  end  result  of  a subcritical  cloud  undergoing  ambipolar 
diffusion  is  an  isolated,  low  mass  star. 

In  the  bimodal  model,  the  magnetic  field  strength  is  the  determining  factor 
between  isolated  and  clustered  star  formation.  However,  observational  studies  have 
failed  to  identify  any  magnetically  subcritical  clouds  (Crutcher  1999).  Furthermore, 
observations  indicate  that  typically  only  a few  percent  of  a cloud’s  mass  is  convert- 
ed into  stars  (Leisawitz  et  al.  1989).  Mechanisms  are  needed  that  can  lend  global 
support  to  a cloud  while  facilitating  regions  of  local  collapse  within  the  cloud.  The 
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dynamic  nature  of  the  ISM  must  be  considered.  The  free-fall  time  for  gravitational 
collapse  is  greatly  modified  if  the  initial  state  of  the  gas  is  not  static.  Convergent 
supersonic  velocity  fields  can  reduce  the  mass  needed  for  collapse  to  a value  below  the 
classical  Jeans  mass  (Hunter  1979,  Hunter  and  Fleck  1982).  Recently,  Klessen  et  al. 
(2000)  performed  simulations  which  demonstrate  that  random,  supersonic  motions 
can  support  a cloud  against  global  collapse,  while  colliding  streams  within  the  cloud 
give  rise  to  regions  unstable  to  local  collapse.  Thus,  star  formation  can  occur  on  small 
scales  within  a cloud  even  though  the  cloud  is  supported  against  global  collapse.  A 
new  model  for  the  ISM  has  emerged  that  further  emphasizes  the  role  of  turbulence. 
This  theory  holds  that  clouds  are  density  enhancements  created  by  the  interaction  of 
larger  scale  turbulent  flows  (Elmegreen  1993b,  Fleck  1996,  Ballesteros-Paredes  et  al. 
1999a).  The  self-similar  nature  of  clouds  is  then  explained.  Convergent  flows  within 
a cloud  give  rise  to  density  enhancements  that  appear  to  form  smaller,  denser  clouds. 
Ballesteros-Paredes  et  al.  (1999b)  speculate  that  clouds  and  stars  may  form  on  time- 
scales  far  shorter  than  anticipated  by  quasi-static  theories.  They  suggest  that  the 
Taurus  molecular  cloud  formed  and  evolved  into  a star  forming  region  in  less  than 
three  million  years. 

1.2  Hydrodynamic  Simulations  of  the  ISM 

The  role  of  turbulence  and  hydrodynamic  instabilities  in  shaping  the  ISM  has 
been  investigated  for  the  last  half  century.  In  one  early  work,  von  Weizsacker  (1951) 
considered  the  evolution  of  stars  and  galaxies  in  a turbulent  and  compressible  medi- 
um. Chandrasekhar  and  Munch  (1952)  studied  the  ISM  as  a “fluctuating  continuous 
distribution  of  matter”.  In  both  of  these  visionary  works,  the  authors  abandoned 
the  idea  that  the  ISM  can  be  characterized  by  discrete  clouds.  They  recognized  the 
importance  of  dynamical  forces  in  the  ISM. 

Due  to  its  highly  nonlinear  nature,  turbulence  cannot  be  modeled  analytically. 
Further  complications  arise  for  astrophysical  systems.  Unlike  terrestrial  laboratory 
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turbulence,  interstellar  turbulence  is  highly  supersonic,  magnetized,  self-gravitating, 
and  contains  embedded  sources  (Elmegreen  1993a).  And  whereas  terrestrial  turbu- 
lence is  driven  by  the  injection  of  energy  only  on  large  scales,  energy  is  injected 
into  interstellar  turbulence  on  many  scales.  Galactic  shear  and  superbubbles  induce 
turbulent  motions  on  large  scales  (>  Ikpc).  Mechanisms  which  operate  at  mod- 
erate scales  (>  lOOpc)  include  supernovae  and  expanding  H II  regions.  On  small 
scales  within  molecular  clouds  (<  10 pc),  stellar  winds  and  protostellar  outflows  im- 
part kinetic  energy  to  the  surrounding  environment.  Numerical  modeling  is  needed 
to  investigate  the  complex  forces  at  work.  Advances  in  computing  power  have  led 
to  an  explosion  in  the  number  of  simulations  of  the  ISM.  (A  recent  review  is  given 
by  Vazquez- Semadeni  et  al.  2000.)  The  first  simulations  of  colliding  clouds  (Stone 
1970a,  1970b)  required  a high  degree  of  symmetry.  Far  more  complicated  simulations 
are  now  possible.  Although  the  formation  of  a comprehensive  model  for  the  ISM  is 
still  infeasible  today,  the  range  of  physical  processes  under  consideration  is  expand- 
ing. Researchers  now  recognize  the  importance  of  thermodynamic  processes  in  the 
ISM.  For  instance,  approximating  the  ISM  as  an  isothermal  gas  is  inappropriate  in 
most  situations  (e.g.  Hunter  et  al.  1986,  Tohline  et  al.  1987,  Vazquez-Semadeni  et 
al.  1996).  The  development  of  magnetohydrodynamics  codes  has  allowed  the  inclu- 
sion of  magnetic  fields  and  plasma  flows.  As  computational  simulations  continue 
to  advance,  astronomers  must  remember  to  carefully  consider  the  physical  processes 
driving  the  dynamics  and  evolution  of  the  ISM.  Fully  understanding  the  physics  of 
the  ISM  requires  careful  isolation  and  study  of  each  process  in  turn.  Only  then  can 
processes  be  coupled  together  in  efforts  to  form  models  which  are  more  in  line  with 
the  observational  picture  of  the  ISM. 

The  ISM  is  susceptible  to  a wide  range  of  hydrodynamic  instabilities.  Insta- 
bilities fall  largely  into  two  categories,  global  and  interfacial.  Global  instabilities  act 
throughout  an  entire  volume  of  gas.  Examples  include  the  Jeans  instability  (Jeans 
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1929),  thermal  instability  (Weymann  1960,  Field  1965),  and  Parker  instability  (Park- 
er 1966).  Interfacial  instabilities  act  locally  at  an  interface.  Examples  include  the 
Rayleigh- Taylor  instability,  Kelvin-Helmholtz  instability,  and  thermal-dynamic  in- 
stability (Hunter  and  Whitaker  1989).  The  effort  to  understand  how  instabilities 
develop  is  complicated  when  conditions  give  rise  to  multiple  instabilities  at  the  same 
time.  Recreating  astrophysical  conditions  on  Earth  is  very  difficult  at  best.  Numer- 
ical simulations  are  the  best  tools  for  isolating  and  studying  instabilities.  Extensive 
numerical  studies  have  been  performed  to  examine  the  Rayleigh-Taylor  instability 
(e.g.  Youngs  1984,  Jun  et  al.  1995)  and  the  Kelvin-Helmholtz  instability  (e.g.  Mur- 
ray et  al.  1993,  Frank  et  al.  1996,  Vietri  et  al.  1997).  However,  some  instabilities 
have  not  received  a great  deal  of  attention.  Two  such  instabilities  provide  the  focus 
for  this  dissertation.  The  first  is  a recently  discovered  (Hunter  et  al.  1997)  insta- 
bility which  acts  across  an  interface  of  discontinuous  density  and  is  driven  solely  by 
self-gravity.  The  second  is  a dynamic  instability  which  occurs  when  supersonic  gas 
streams  collide.  This  instability  is  driven  by  ram  pressure  and  strong  cooling  (Hunter 
et  al.  1986,  Stevens  et  al.  1992).  At  issue  is  whether  interacting  streams  give  rise  to 
collapsing  regions.  If  so,  then  one  epoch  of  star  formation  may  trigger  the  subsequent 
birth  of  additional  stars  (Elmegreen  and  Lada  1977).  Both  instabilities  are  described 
in  great  detail  in  later  chapters. 

1.3  Thesis  Structure 

This  dissertation  is  divided  into  six  chapters.  The  first  is  this  introductory 
chapter.  In  chapter  2,  I describe  the  numerical  methods.  In  chapter  3,  I present  the 
first  computational  study  to  focus  on  interfacial  instabilities  driven  by  self  gravity. 
Chapters  4 and  5 are  devoted  to  the  study  of  colliding  gas  streams  and  the  instabilities 
that  develop.  In  chapter  4,  I concentrate  on  models  in  Cartesian  geometry  in  which 
the  role  of  molecular  cooling  is  emphasized.  In  chapter  4,  I discuss  cylindrical  models 
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with  a focus  on  the  formation  of  collapsing  clumps.  I conclude  this  work  in  chapter 
6 with  a summary  of  the  results  and  some  remarks  about  future  considerations. 


CHAPTER  2 

COMPUTATIONAL  TECHNIQUES 


The  study  of  nonlinear  phenomena  requires  the  use  of  numerical  modeling. 
For  problems  concerning  hydrodynamic  instabilities  in  the  ISM,  the  equations  gov- 
erning the  flow  of  gas  are  discretized  and  solved  on  a computer.  The  computer  code 
advances  the  hydrodynamic  state  of  the  gas  in  time.  The  code  used  for  this  study  is 
the  Computational  Fluid  Dynamics  LIBrary  (CFDLIB)  created  at  the  Los  Alamos 
National  Laboratory  (LANL).  In  order  to  model  the  environment  of  interstellar  space, 
specialized  subroutines  are  needed  to  simulate  cooling  processes  and  self-gravitational 
forces.  These  routines  were  written  and  implemented  to  run  within  the  structure  of 
CFDLIB.  Figure  2.1  shows  the  main  subroutines  as  called  by  CFDLIB  along  with  a 
brief  description  of  each.  An  input  file  is  used  to  specify  the  computational  param- 
eters, grid  parameters,  and  initial  values  for  the  state  variables.  CFDLIB  generates 
a grid  and  begins  the  computational  cycle.  Each  cycle  consists  of  several  distinct 
phases.  Subsequent  sections  of  this  chapter  describe  the  hydrodynamic,  self-gravity, 
and  cooling  portions  of  the  code. 

2.1  The  Hydrodynamic  Code 

Hydrodynamic  calculations  are  performed  using  CFDLIB,  a package  of  FOR- 
TRAN codes  from  LANL  capable  of  solving  a wide  range  of  computational  fluid 
dynamics  problems.  CFDLIB  is  well  suited  for  handling  flows  of  any  speed,  from 
low  velocity,  incompressible  flows  to  high  velocity,  compressible  flows.  The  programs 
within  CFDLIB  are  connected  by  a common  data  structure  and  numerical  method, 
which  allows  for  easy  modifications  and  additions  to  the  code  library.  CFDLIB  uses 
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begin  CFDLIB 


rdinput 


Sets  default  values  for  all  quantities  and 
reads  input  file 


O 


Repeat 
cycle 
until 
end  time 

reached 

or 

time  step 
becomes 
less  than 
a specified 
value 


Writes  output  data  at  specified  times 


Drives  the  primary  and  Lagrangian  phases 
of  the  calculation  cycle 


Modifies  temperature  according  to  cooling  model 


Remaps  solution  to  stationary  grid- 
the  Eulerian  phase  of  the  cycle 

Solves  Poisson’s  eq^uation  for  the  gravitational 
potential  and  calculates  gradients  of  the  potential 


Updates  state  variables  & geometric  arrays 
and  finds  new  time  step  for  next  cycle 


Writes  final  output  and 
timing  information 


Figure  2.1:  Chart  of  subroutines  called  by  CFDLIB  main  program. 
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Finite- Volume  computational  schemes  on  structured  grids  in  which  all  state  variables 
are  cell-centered.  In  Finite- Volume  methods,  the  integral  form  of  the  conservation 
equations  are  solved  over  a grid  in  which  each  cell  is  treated  as  a control  volume 
(Kashiwa  et  al.  1994).  Fluid  flow  is  modeled  by  tracking  the  fluxes  of  material,  mo- 
mentum, and  energy  between  adjacent  cells.  Features  of  CFDLIB  allow  the  user  to 
choose  boundary  conditions,  geometry,  an  equation  of  state,  and  whether  to  include 
magnetic  flelds  in  the  calculation.  Calculations  may  be  done  in  two  dimensions  or 
in  three  dimensions,  although  run  parameters  are  more  limited  in  three  dimensions. 
Also,  CFDLIB  is  capable  of  performing  multi-fluid  calculations.  One  of  the  advan- 
tages of  CFDLIB  over  other  codes  is  the  incorporation  of  the  full  viscous  stress  tensor 
into  the  equations  of  hydrodynamics.  Other  codes  rely  on  an  artiflcial  viscosity  term 
to  represent  distortions  of  the  scalar  pressure.  The  full  tensor  gives  a more  realistic 
solution,  especially  in  shocked  regions. 

The  underlying  method  used  by  CFDLIB  is  called  split  Arbitrary-Lagrangian- 
Eulerian  (ALE).  This  method  allows  the  grid  to  move  in  an  arbitrary  way.  The  user 
may  choose  to  use  an  Eulerian  mesh  of  spatially  fixed  grid  points,  a Lagrangian 
mesh  that  moves  with  the  fluid,  or  a combination  of  the  two.  Regardless  of  the 
choice,  each  computational  cycle  is  split  into  three  phases.  In  the  primary  phase,  the 
equation  of  state  is  used  to  solve  for  the  pressure,  and  auxiliary  variables  that  are  not 
explicitly  time  dependent  are  calculated.  Next,  the  effects  of  the  physical  processes 
are  computed  in  the  Lagrangian  phase.  Lastly,  in  the  Eulerian  phase,  the  state 
variables  are  advected  and  mapped  to  a common  control  volume.  CFDLIB  combines 
the  ALE  method  with  a semi-implicit  scheme  known  as  ICE  (Implicit  Continuous- 
fluid  Eulerian).  The  ICE  scheme  is  semi-implicit  in  that  the  implicitness  is  restricted 
to  the  Lagrangian  phase  of  the  calculation.  The  ALE  method  was  developed  for 
modeling  compressible  fluids,  but  the  addition  of  the  ICE  technique  allows  the  code 
to  work  for  any  flow  speed  (Hirt  et  al.  1974).  In  other  words,  CFDLIB  is  stable  for 
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any  value  of  the  Courant  number  based  on  the  sound  speed.  Following  the  description 
of  Addessio  et  al.  (1992),  I present  here  the  basic  ideas  of  the  finite  volume  method 
at  the  heart  of  CFDLIB.  Throughout  this  chapter,  vector  quantities  are  denoted  by 
boldface  characters. 

Neglecting  magnetic  fields,  the  equations  of  fluid  dynamics  which  describe 
the  conservation  of  mass,  momentum,  and  energy  can  be  expressed  in  the  familiar 
differential  form 

^ + V-  (pu)  = 0,  (2.1) 

^ + V • (puu)  = - Vp  + V • n + pF,  (2.2) 

^ + V • (pEu)  = -V-pu  + V-  n-  u + pF-u-V-h  + (5.  (2.3) 

In  these  equations,  p is  the  mass  density,  u the  fluid  velocity,  and  p the  fluid  pressure. 

The  total  energy  per  unit  mass  E is  given  hy  E = e + (l/2)u-u,  where  e is  the  internal 
energy  per  unit  mass.  The  pressure  is  determined  from  p and  E via  the  equation  of 
state.  The  stress  tensor,  II,  describes  any  viscosity  or  turbulence.  For  some  problems, 
additional  physics  may  be  needed  as  described  by  a body  force  per  unit  mass  F,  a 
heat-flux  vector  h,  or  an  energy  release  rate  per  unit  volume  Q.  Gravitational  forces 
are  included  in  F.  Equation  (2.2)  is  known  as  the  Navier-Stokes  equation  for  viscous 
flows.  In  the  Finite- Volume  method,  the  integral  form  of  the  equations  is  used.  The 
integral  form  expresses  the  idea  that  a conserved  quantity  (mass,  momentum,  or 
energy)  may  only  change  within  a particular  control  volume  by  an  amount  equal 
to  the  net  flux  of  that  quantity  across  the  bounding  surface  of  the  control  volume. 
The  values  within  a cell  are  updated  by  keeping  track  of  the  amount  of  material, 
momentum,  and  energy  entering  and  leaving  the  cell.  When  expressed  in  integral 
form,  the  equations  are 


d 

dt 


pdV  + 


p(u  - u,„)  • ndS  = 0, 


(2.4) 
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d 

dt 


/ pudV  + / pu(u  — Urn)  • ndS  = — ; 

J Srr.it)  J Srr.it) 


pndS 


+ 


[ n-UdS+  [ pFdV, 

JSrr.it)  JVrr.it) 


(2.5) 


£ 

dt 


/ pEdV  + / pE{u  — u^)  • ndS  = — pu  • nc?5  + 

yK„(t)  ys,n(t) 


[ n U udS-  [ h ndS+  [ pF  ■ udV  [ QdV.  (2.6) 

JSrr.it)  JSrr.it)  J Vrr.it)  J Vrr.it) 

Each  cell  acts  as  a control  volume  Vm{t)  with  a surface  of  Sm{t)  and  outward  unit 
normal  n.  The  control  volume  moves  with  velocity  u^-  The  volume  in  each  cell  is 
conserved,  so  that 

dVm{t) 


dt 


= j ^ 

JSrr.it) 


.ndS. 


(2.7) 


In  the  Lagrangian  limit,  the  control  volume  moves  at  the  same  velocity  as  the  fluid  so 
u — Urn  = 0.  This  simplifles  the  equations.  For  example,  the  equation  of  continuity, 
equation  (2.4),  becomes 

= (2.8) 


The  Lagrangian  control  volume  is  designated  V^.  The  ICE- ALE  method  takes  ad- 
vantage of  the  simplifled  equations  by  updating  the  variables  in  the  Lagrangian  limit 
and  then  remapping  the  solution  onto  the  control  grid. 

CFDLIB  implements  the  ICE-ALE  method  in  a way  that  is  modified  from  the 
original  algorithm  of  Hirt  et.  al  (1974).  The  use  of  the  Lagrangian  and  Eulerian  phases 
remains  the  same,  but  a new  primary  phase  of  the  calculation  is  added  to  improve 
the  efficiency  of  the  code  and  generalize  the  code  for  multiphase  flows.  Following 
Kashiwa  et  al.  (1994),  I summarize  the  CFDLIB  implementation  of  the  ICE-ALE 
method  beginning  with  the  Eulerian  phase  and  proceeding  backwards  to  the  primary 
phase.  This  way  the  motivation  for  the  primary  phase  shall  be  made  clear.  The  main 
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variables  in  CFDLIB  are  the  meiss  density  p,  the  velocity  u,  and  the  temperature 
T.  The  choices  of  velocity  over  momentum  and  temperature  over  energy  help  to 
avoid  generating  numerical  noise.  The  code  need  not  divide  by  the  total  mass  to 
define  the  velocity  nor  convert  from  total  energy  to  temperature.  It  is  especially 
important  to  resolve  the  temperature  as  accurately  as  possible  for  problems  with 
chemical  reactions.  In  an  ideal,  polytropic  gas,  the  temperature  and  specific  internal 
energy  are  related  simply  as  e = c„T,  with  c„  defining  the  specific  heat  at  constant 
volume.  The  equation  of  state  relates  the  temperature,  density,  and  pressure  so  that 
any  two  can  be  used  to  determine  the  third.  For  an  ideal  gas,  the  equation  of  state  is 

P = —PT,  (2.9) 

where  7 is  the  ratio  of  specific  heats,  fj,  the  molecular  weight,  and  R the  gas  constant, 
(i?  = 8.314 X 10^  mol~^).  Both  the  ideal  gas  law  and  the  polytropic  relation, 

p ca  , are  used  throughout  this  study. 

Although  CFDLIB  has  been  generalized  to  handle  multiple  species  of  fiuid,  I 
describe  the  single  fluid  version  to  simplify  the  notation.  The  variable  q is  used  to 
represent  any  state  variable.  In  this  description,  q is  considered  the  average  value 
within  a control  volume  V so  that  fy  qdV  = qV.  The  advection  equation  that  relates 
the  total  change  in  g on  a computational  mesh  moving  at  to  the  change  in  q due 
to  the  physical  processes  within  V is 

+ J n • g(u  - Um)dS  = (2.10) 

The  right  hand  side  represents  the  change  due  to  physical  processes,  which  is  the 
solution  from  the  Lagrangian  phase  of  the  calculation.  Using  superscripts  to  denote 
the  time  step,  equation  (2.10)  is  discretized  and  written 


qU+lyn+l  _ ^nyn  ^ AtA(g)  = q^V^  ~ 


(2.11) 
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The  terms  cancel  leaving  a relationship  between  the  updated  value  on  the  mesh 
and  the  Lagrangian  value.  The  advection  operator  AtA{q)  is  defined  as 

AtA{q)  = ^ AtSiiu*  - u^),((g))^  (2.12) 

i 

The  sum  is  taken  over  all  surfaces  defining  y",  with  Si  representing  the  area  times 
outward  normal  of  each  side.  The  face-centered  fluxing  velocity  u*  is  one  of  the  main 
results  of  the  primary  phase.  The  bracketed  term  is  an  upwind-centered  value  of  q, 
expanded  in  a spatial  Taylor  series  to  a point  just  upstream  of  i.  The  Taylor  series 
is  second-order  in  space  and  uses  gradients  that  are  limited  using  the  method  of  van 
Leer  (1977).  The  Lagrangian  volume  and  the  new  mesh  volume  1^"+^  are  both 
updated  to  reflect  changes  to  the  velocity  fields, 

+ At^(S-u*)i  (2.13) 

i 

yn+i  = v^  + At  ^(S  • u™)i  (2.14) 

i 

The  updated  volumes  are  used  along  with  equation  (2.11)  to  arrive  at  the  final  equa- 
tions for  the  updated  state  variables, 


p"+i  = [p^V^  - At^(p")]/F"+^ 

(2.15) 

^ ^^LyL  _ Atyl(u")]/F"+^ 

(2.16) 

rpn+l  ^ ^j^LyL  _ Afyl(r")]/F"+\ 

(2.17) 

The  values  of  p^,  u^,  and  are  determined  in  the  Lagrangian  phase. 

In  the  Lagrangian  phase,  the  effects  of  the  physical  processes  are  calculated. 
The  Lagrangian  values  for  density,  velocity,  and  temperature  are  determined  in  a 
manner  similar  to  that  of  equations  (2.15-2.17),  namely 

= [p"F"  + Am]/F^  (2.18) 

= u"[l  — (Am/m'^j  + {Amu)/m^ 


(2.19) 
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TL  ^ _ (en/c„T")(Am)/m^]  + (Ame)/m^.  (2.20) 

The  value  of  Am  reflects  the  exchange  of  mass  between  different  species  of  fluid.  For 
the  case  of  a single  fluid,  the  Am  term  vanishes,  leaving  a simplifled  set  of  equations. 


= P^V'^IV^ 

(2.21) 

: u"  -h  (Amu)/m^ 

(2.22) 

T"  + (Ame)/m^. 

(2.23) 

For  the  remainder  of  this  section,  I adopt  the  notation  of  Kashiwa  et  al.  (1994)  and 
use  to  denote  the  discrete  form  of  a cell-centered  integral  over  the  sides  of  a 

computational  cell  which  acts  on  face-centered  quantities.  The  Amu  and  Ame  terms 
contain  the  physics  for  the  problem, 

(Amu) /m^  = (F/m) At  - (p“^ Atl/" V V) /V^  + ( AtV” V"-r) /m^,  (2.24) 

{Ame)/m^  = {pc)~^{p^Ap)/V^.  (2.25) 

The  quantity  c is  the  sound  speed  (c^  = 'yP/p).  These  equations  do  not  include 
turbulence  and  are  valid  only  for  single-phase  fluids.  The  momentum  equation,  e- 
quation  (2.24),  contains  the  momentum  change  due  to  body  forces,  pressure  forces, 
and  forces  from  stress.  Gravitational  forces  are  included  in  the  body  force  term  F, 
and  T represents  the  stress.  The  energy  equation,  equation  (2.25),  looks  deceptively 
simple.  However,  the  expression  for  Ap  is  rather  complicated.  The  quantities  p*  and 
Ap  are  determined  in  the  primary  phase. 

The  primary  phase  of  the  calculation  is  responsible  for  calculating  quantities 
which  are  not  explicit  functions  of  time.  The  equation  of  state,  equation  (2.9)  for 
an  ideal  gas,  is  evaluated  to  find  the  pressure.  Also  calculated  in  the  primary  phase 
are  the  auxiliary  quantities  used  in  the  Lagrangian  and  Eulerian  phases.  The  re- 
quired quantities  are  the  face-centered  fluxing  velocity  u*  (used  in  equations  2.12 
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and  2.13),  the  face-centered  equilibration  pressure  p*  (used  in  equation  2.24),  and 
the  cell-centered  change  in  the  equilibration  pressure  Ap  (used  in  equation  2.25).  For 
ease  of  notation,  let  V-^  represent  a face-centered  gradient  operating  on  surrounding 
cell-centered  data.  Also,  let  v represent  the  specific  volume  (volume  per  unit  mass). 
In  models  involving  multiple  fluids,  the  value  of  v must  be  tracked  for  each  individual 
fluid.  In  the  single  fluid  case,  v is  simply  the  reciprocal  of  the  density,  v = For 
further  convenience,  the  quantity  is  defined  as  = l/{pi  + Pr),  with  the  sub- 
scripts I and  r denoting  the  values  in  the  cell  to  the  left  or  right  of  the  cell  face  being 
considered. 

The  equilibration  pressure  is  defined  for  multiple  fluids  as  the  pressure  that 
enables  specified  masses  of  separated  materials  to  fill  an  entire  volume  with  no  ongoing 
compression  or  expansion.  In  the  single  fluid  limit,  the  equilibration  pressure  is  simply 
the  pressure  at  a cell  face,  as  opposed  to  the  cell-centered  pressure  calculated  from 
the  equation  of  state.  The  change  in  the  equilibration  pressure  is  determined  from 

[V^{v/c^)  - • vN^]Ap  = • u,  (2.26) 


where 

u=(  + ^ g^^/2. 

V Pi  + Pt  J 


(2.27) 


The  vector  g represents  the  acceleration  due  to  all  body  forces,  not  just  gravity.  Find- 
ing Ap  means  solving  a linear  system  of  equations  that  involves  every  computational 
cell.  Once  this  has  been  done,  the  fluxing  velocity  is  defined  as 


u*  = I ) _ t;/AtVV  + gAt/2.  (2.28) 

V Pi  + Pr  J 

Notice  that  the  equation  for  the  fluxing  velocity  is  the  same  as  equation  (2.27)  except 
that  p”  has  been  replaced  by  p^  = p^  + Ap.  This  is  one  of  the  major  factors  in  the 
stability  of  the  ICE  method.  Linear  stability  is  assured  when  p^  is  used  in  the 
equation  for  the  pressure  because  p^  provides  a good  estimate  of  the  advanced-time 
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pressure.  The  equilibration  pressure  p*  used  in  the  equation  of  motion  is  given  by 
. _ f VrPr  + i^Pr  ~ ^ [^r(^r  ~ U*)  - Ui{u*  - Uj)] 

\ Vr  +Vi  ) At  Vr  + U;  Ur  + '^^1 

A flux  limiter  0 is  applied  to  the  second  two  terms  in  order  to  minimize  diffusion. 
The  idea  is  to  limit  any  non-physical  spatial  distribution  in  the  velocity  field  that 
arises  out  of  numerical  noise.  The  limiter  has  a value  between  zero  and  one,  with  the 
value  tending  towards  zero  when  the  velocity  held  is  smooth. 

The  overall  algorithm  of  the  ICE- ALE  algorithm  as  implemented  by  CFDLIB 
begins  by  calculating  the  equation  of  state  and  auxiliary  variables  in  the  primary 
phase.  Next,  the  physical  processes  are  computed  in  the  Lagrangian  phase.  Finally, 
the  solution  is  advected  onto  the  control  mesh  in  the  Eulerian  phase.  If  the  user 
chooses  to  run  the  code  as  purely  Lagrangian,  the  Eulerian  step  is  unnecessary.  The 
time  interval  between  cycles  is  limited  by  the  Courant  timestep  as  determined  from 
the  greater  of  the  flow  speed  or  the  sound  speed.  Examples  of  model  runs  done  with 
CFDLIB  are  given  by  Kashiwa  et  al.  (1994).  What  follows  is  an  example  of  a purely 
hydrodynamic  model  done  with  CFDLIB  that  is  applicable  to  astrophysical  scales. 

When  a fluid  is  placed  on  top  of  a more  tenuous  fluid  in  a gravitational  field 
(such  that  the  denser  fluid  is  accelerated  into  the  tenuous  fluid),  the  interface  exists 
in  an  unstable  equilibrium.  If  the  interface  is  disturbed,  the  denser  fluid  will  displace 
the  lighter  fluid,  ultimately  resulting  in  the  inversion  of  the  original  configuration. 
The  name  given  to  this  phenomenon  is  the  Rayleigh- Taylor  (RT)  instability.  For  an 
incompressible  gas,  the  linear  growth  of  an  RT  instability  is  well  documented  (e.g. 
Chandrasekhar  1961).  For  a perturbation  proportional  to  the  growth  rate 

is  given  by 


where  uj  is  the  growth  rate,  k the  horizontal  wavenumber,  and  g the  gravitational 


(2.30) 


acceleration  in  the  vertical  direction.  Figure  (2.2)  shows  the  configuration  of  the  RT 
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Figure  2.2:  Configuration  of  the  Rayleigh-Taylor  simulation. 

test  problem.  The  simulation  is  performed  on  an  80x80  Cartesian  grid  with  a fixed 
cell  size  of  Ax  = Ay  = 1.25  x 10^®  cm.  The  grid  is  constrained  by  reflective  boundary 
conditions  on  all  four  sides.  The  grid  is  treated  as  though  filled  with  molecular 
hydrogen  gas  with  molecular  weight  = 2.0  and  a ratio  of  specific  heats  7 = 1.4. 
The  number  density  and  temperature  below  the  interface  are  rii  = 100  cm“^  and 
Ti  = 20  K,  respectively.  The  number  density  and  temperature  above  the  interface  are 
ri2  = 200  cm~^  and  T2  — 10  K,  respectively.  The  pressure  is  continuous  at  the  interface. 
The  temperature  remains  constant  throughout  each  region,  but  the  density  falls  off 
as 

Iv-yWI 

n{y)i,2  = «(0)i,2e  "1-2  , (2.31) 

where  y(0)  = 5.0  x 10^^  cm  is  the  position  of  the  interface.  The  gravitational  scale 
height  is 

H = (2.32) 

IJLmng 

where  ke  = 1.38067  x 10'^®  cry is  Boltzmann’s  constant,  and  ttih  — 1.6735  x 
10“^^  g is  the  mass  of  a hydrogen  atom.  The  gravitational  acceleration  is  chosen 
so  that  the  density  falls  to  0.25e^/^n(0)  over  region  2,  resulting  in  g = —1.03125  x 


Region  2 


T=  10  K 
n=200  ciri’ 


Region  1 

T=20K 
n=100  cni’ 


10  cm  = 0.32  pc 
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10  ^°cms  At  t = 0 a perturbation  is  applied  to  the  vertical  velocity  at  the  interface 
such  that 

2ttx 

where  C2  is  the  sound  speed  in  region  2,  = 'yksT2/ A small  fraction  of 

the  sound  speed  is  chosen  as  the  maximum  perturbation  amplitude  in  order  to 
approximate  an  incompressible  fluid.  The  chosen  wavelength  is  one- fourth  the  to- 
tal grid  size,  A = 2.5  x 10^^  cm,  which  corresponds  to  a horizontal  wave  number 
k = 27t/A  = 2.51  X 10“^^  cm~h  Solving  for  the  linear  growth  rate  using  equation 
(2.30)  yields  a;  = 2.94  x 10~^'*  Inverting  the  growth  rate  results  in  the  e-folding 
time  for  the  problem,  = = 3.40  x 10^^  s.  Figure  (2.3)  shows  the  density  distribu- 

tion at  4,6,8,  and  10  e-folding  times.  Early  in  the  simulation,  the  visible  growth  of  the 
instability  is  small  due  to  the  large  cell  size,  small  velocity  perturbation,  and  small 
gravitational  acceleration.  The  initial  wavelength  is  preserved  well  into  the  nonlinear 
phase  of  the  instability.  The  boundary  conditions  give  rise  to  some  asymmetry  near 
the  left  and  right  boundaries.  As  the  instability  evolves,  the  downward  velocity  of  the 
heavy  fluid  increases.  At  late  times,  the  velocity  shear  is  great  enough  to  give  rise  to 
mushroom  caps  due  to  the  Kelvin-Helmholtz  instability.  These  results  compare  well 
to  laboratory  scale  calculations  by  Jun  et  al.  (1995).  More  detailed  RT  simulations 
are  presented  in  the  next  chapter. 

2.2  Self-Gravity  Routines 

On  astrophysical  scales,  gravitational  forces  play  an  enormous  role  in  the 
evolution  of  gas  clouds.  In  order  for  a gas  cloud  to  evolve  into  a protostar,  gravity 
must  at  some  point  overwhelm  the  pressure  forces  acting  against  collapse.  To  study 
the  effects  of  self-gravity  upon  interfacial  instabilities,  I built  a gravity  solver  to  be 
implemented  within  the  structure  of  CFDLIB.  This  section  describes  the  methods 
used  for  the  gravity  solver  and  shows  the  test  problem  of  a collapsing  sphere. 


Vy{x,y)  = .001c2sin 


— lu- 
e A 


v(o)| 


(2.33) 
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Figure  2.3:  Rayleigh- Taylor  instability  growth.  Contours  show  the  number  density 
of  molecular  hydrogen  at  a)4te;  b)6te;  c)8te]  and  d)10te-  The  e-folding  time  = 
3.40  X 10^^  s. 
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Forces  arising  from  self-gravity  are  easily  added  into  the  ICE- ALE  procedure 
by  incorporating  them  into  the  body  force  terms  in  the  primary  and  Lagrangian  steps. 
The  forces  themselves  are  given  by  the  gradient  of  the  self-gravitational  potential, 
f = V$.  It  is  simpler  and  quicker  to  solve  for  the  scalar  potential  than  the  vector 
forces.  The  self- gravitational  potential  $ for  a fluid  is  found  by  solving  the  Poisson 
equation 

= -AttGp,  (2.34) 

where  G is  Newton’s  Gravitational  constant  (G  = 6.6742  x cm^g~^s~^).  Where- 
as the  equations  of  fluid  dynamics,  equations  (2. 1-2.3),  are  hyperbolic  equations,  the 
Poisson  equation  is  an  elliptic  equation.  Elliptic  equations  are  also  called  boundary 
value  problems  (both  hyperbolic  and  parabolic  equations  are  initial  value  problems) 
and  require  different  methods  to  solve  than  hyperbolic  equations.  Poisson’s  equation 
contains  no  time  dependence;  the  potential  value  in  each  cell  throughout  the  mesh 
must  be  solved  for  simultaneously.  Finding  the  solution  comes  down  to  solving  a set 
of  nonlinear  algebraic  equations  by  linearization  and  iteration.  The  solution  is  very 
sensitive  to  the  boundary  conditions.  Many  methods  exist  to  solve  elliptic  equations, 
such  as  Fourier  techniques,  conjugate  gradient  methods,  or  multigrid  methods.  The 
solver  implemented  in  CFDLIB  is  a multigrid  solver.  The  multigrid  method  is  as  fast 
as  other  methods,  and  it  retains  more  generality  in  choosing  boundary  conditions 
and  cell  sizes.  The  name  of  the  employed  multigrid  code  is  MUDPACK  (Adams 
1989,1991),  a freely  available  program  developed  by  John  Adams  at  the  National 
Center  for  Atmospheric  Research.  MUDPACK  is  a collection  of  Fortran  subroutines 
created  to  make  multigrid  methods  available  in  a user  friendly  form.  The  self-gravity 
routine  added  to  CFDLIB  consists  of  parts  of  MUDPACK  along  with  code  to  allow 
communication  between  CFDLIB  and  MUDPACK.  The  potential  values  are  calcu- 
lated after  the  state  variables  have  been  remapped  to  the  Eulerian  grid  to  ensure  that 
the  cell  size  is  a constant  in  each  direction. 
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Multigrid  methods  combine  the  use  of  basic  iteration  techniques  with  a process 
of  stepping  the  solution  down  to  coarser  grids  to  speed  convergence.  With  iteration 
techniques,  a sequence  of  approximate  solutions  is  generated  until  convergence  to 
a solution  is  achieved.  As  an  example,  consider  the  one-dimensional,  second  order 
discrete  form  of  the  Poisson  equation, 

$i_i  — 2$j  -t“  = — AttG pih^ , (2.35) 


where  h is  the  cell  size.  Letting  bj  = —AirCpib?,  this  equation  can  be  written  in  matrix 
form  as  Ax  = b with  the  definitions. 
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The  matrix  A is  symmetric,  tridiagonal,  and  sparse.  Iterative  techniques  exploit  the 
sparse  nature  of  the  matrix.  The  simplest  iterative  method  is  the  Jacobi  method. 
This  method  starts  off  with  an  initial  guess  for  x.  All  the  components  of  x are  then 
updated  simultaneously  using  the  simple  relation, 

= (5<  - ^ aijxf  Vai.i-  (2-36) 

j¥=i 


Since  A is  sparse,  most  of  the  aij  terms  are  zero,  so  each  value  of  Xi  is  updated  as 
a function  only  of  its  closest  neighbors.  The  iteration  continues  until  convergence  is 
achieved.  Notice  that  in  the  Jacobi  method,  each  updated  value  is  a function 

only  of  values  from  the  previous  iteration.  The  order  in  which  the  updated  values 
are  calculated  does  not  matter.  The  Gauss-Seidel  method  is  similar  to  the  Jacobi 
method,  except  that  updated  values  are  used  as  soon  as  they  are  calculated.  The 
order  in  which  the  components  of  x are  updated  becomes  important.  For  example, 
values  can  be  updated  according  to  the  relation 


J<i 


j>i 


(2.37) 


23 


The  Gauss-Seidel  method  requires  less  storage  than  the  Jacobi  method,  and  will 
converge  with  fewer  iterations  for  most  problems.  The  direction  of  the  iteration 
need  not  be  the  same  each  time.  In  the  symmetric  Gauss-Seidel  method  the  up- 
date order  is  alternated  between  the  ascending  and  descending  directions.  Another 
variation  involves  updating  the  even  and  odd  components  in  alternating  order.  For 
the  tridiagonal  matrix  A,  the  even  components  are  updated  as  a function  of  the 
odd  components,  and  vice  versa.  This  process  is  called  the  red-black  Gauss-Seidel 
method.  More  advanced  iterative  methods  such  as  successive  over- relaxation  and 
weighted  function  methods  exist  as  well,  but  with  multigrid  techniques  the  simplest 
iterative  techniques  are  sufficient. 

MUDPAGK  uses  the  red-black  Gauss-Seidel  method  as  the  iteration  tech- 
nique. Like  most  iterative  techniques,  the  Gauss-Seidel  method  suffers  from  some 
major  drawbacks.  Gonvergence  is  often  very  slow.  If  the  initial  guess  is  very  wrong, 
convergence  may  never  happen.  Multigrid  techniques  provide  a way  to  overcome 
these  problems.  Iterative  techniques  have  a smoothing  property  (Briggs  1987)  in 
that  they  quickly  eliminate  oscillatory  (high  frequency)  modes  while  leaving  smooth 
(low  frequency)  modes  intact.  Just  a few  iterations  are  sufficient  to  remove  the  os- 
cillatory modes,  but  to  remove  the  smooth  modes  requires  many  more.  The  key  to 
multigrid  techniques  is  that  a smooth  mode  on  one  grid  appears  more  oscillatory  on 
a coarser  grid.  Figure  (2.4)  shows  an  example  of  this  phenomenon.  A sine  wave  that 
is  smooth  on  a 12  point  grid  becomes  more  oscillatory  if  every  other  data  point  is 
transferred  onto  a grid  with  one-half  the  resolution.  In  a multigrid  process,  a problem 
might  begin  on  the  finest  grid.  After  a few  iterations  on  the  fine  grid,  a method  of 
restriction  is  used  to  transfer  the  problem  to  a coarser  (usually  by  a factor  of  2 in 
each  direction)  grid  where  another  sequence  of  iterations  is  performed.  A process  of 
interpolation  can  be  used  to  transfer  the  solution  back  onto  the  finer  mesh.  Typically 
more  than  two  grid  sizes  are  used.  The  problem  is  simply  transferred  from  grid  to 


24 


Figure  2.4:  An  example  of  switching  grid  resolutions  to  alter  the  nature  of  a function. 
A smooth  function  on  a fine  12  point  grid  appears  more  oscillatory  on  a coarser  6 
point  grid. 

grid  as  needed  to  achieve  convergence  of  the  overall  problem  in  a short  time.  In  cases 
in  which  there  is  no  initial  guess  available  to  begin  the  iteration  on  the  finest  grid, 
the  process  can  begin  on  the  coarsest  grid  and  advance  to  finer  and  finer  grids.  A 
starting  initial  value  is  determined  for  each  grid  in  turn.  The  ordering  of  the  grids 
is  shown  in  figure  (2.5)  for  this  procedure.  A sequence  of  cycles  is  used  to  build  up 
to  an  initial  approximation  to  the  solution  on  the  finest  grid.  A full  cycle  down  to 
the  coarsest  grid  and  back  up  to  the  finest  grid  is  then  completed  to  improve  the 
solution.  This  process  is  known  as  the  full  multigrid  V-cycle  (FMV)  and  is  one  of 
the  processes  used  by  MUDPACK.  Other  variations  make  use  of  even  more  cycles  to 
improve  the  solution  on  the  finest  grid.  Because  only  a few  iterations  are  needed  on 
each  grid  the  overall  process  remains  very  efficient. 

MUDPACK  accepts  three  types  of  boundary  conditions:  periodic,  specified 
value  (Dirichlet),  and  mixed  derivative  (Neumann).  The  appropriate  type  of  bound- 
ary condition  is  determined  by  the  nature  of  the  simulation.  Periodic  boundary 
conditions  approximate  an  infinitely  extended  grid.  For  some  problems  periodic 
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Figure  2.5:  A schematic  look  at  the  grid  schedule  in  the  FMV  method.  The  process 
begins  on  the  coarsest  grid  and  ends  on  the  finest  grid.  The  characteristic  cell  size 
of  the  finest  grid  is  h,  and  each  subsequent  step  to  a coarser  grid  results  in  twice  the 
cell  size.  A sequence  of  iterations  is  performed  on  each  grid. 

boundaries  work  very  well,  provided  that  the  grid  is  made  large  enough  that  the 
area  of  interest  is  well  removed  from  the  boundaries.  For  models  that  use  cylindrical 
symmetry,  periodic  boundaries  are  not  applicable  at  the  symmetry  axis  or  outer  radi- 
al boundary.  Therefore,  a model  with  cylindrical  symmetry,  or  Cartesian  symmetry 
if  not  well  approximated  by  an  infinitely  extended  grid,  requires  specified  value  or 
derivative  boundaries.  The  choice  between  specified  value  or  derivative  boundaries  is 
usually  determined  by  practical  considerations.  The  derivative  boundary  conditions 
are  applied  on  all  grids,  whereas  specified  value  boundaries  are  only  specified  on  the 
finest  grid.  When  the  potential  remains  relatively  unchanged  for  long  intervals  of 
time,  an  approximate  analytical  formula  for  the  boundaries  can  be  applied.  This  is 
typically  the  case  when  a simulation  is  started  from  a state  of  hydrodynamic  equilib- 
rium and  the  total  mass  over  the  grid  remains  constant.  A slowly  changing  potential 
enables  the  efficient  use  of  both  specified  value  and  specified  derivative  boundary 
conditions.  On  the  other  hand,  if  the  problem  under  consideration  involves  a chang- 
ing total  mass  or  the  gas  is  highly  compressible,  the  potential  will  vary  greatly  over 
time.  In  this  case,  boundary  values  must  be  approximated  computationally  using 
information  within  the  grid,  and  specified  value  boundaries  are  more  practical  since 
they  need  be  calculated  only  once  per  time  step.  For  two-dimensional  cases  with 
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imposed  cylindrical  symmetry,  the  potential  can  be  specified  for  a boundary  cell  by 
explicitly  summing  the  contributions  from  all  other  cells.  Each  cell  defines  an  annulus 
of  constant  density.  The  potential  in  each  boundary  cell  is  determined  by  breaking 
up  each  annulus  into  small  pieces  and  treating  each  piece  as  a point  mass.  The  sum 
of  potentials  from  all  pieces  is  the  potential  for  the  boundary  cell.  Once  the  boundary 
values  are  specified,  MUDPACK  very  quickly  solves  for  the  potential  over  the  interior 
of  the  grid.  Experiments  have  shown  that  increasing  the  number  of  divisions  for  each 
annulus  beyond  24  results  in  only  a minor  increase  in  accuracy. 

A simple  example  for  which  an  analytic  expression  for  the  potential  exists  is 
the  problem  of  a collapsing  sphere.  A sphere  with  internal  density  p and  temperature 
T will  collapse  if  its  radius  exceeds  the  Jeans  length 


— 


pGp  ' 


(2.38) 


For  a sphere  with  radius  a,  the  interior  potential  at  a distance  r from  the  center  is 


4>(r)  (2/3)7rGp(3a^  — r^). 


(2.39) 


The  external  potential  for  r > a is 

4>(r)  = (4/3)7rGpa^r“h 


(2.40) 


To  calculate  the  potential  numerically,  a 49x97  grid  is  created  with  cylindrical  geom- 
etry. The  cell  size  is  5 x 10^^  cm  for  both  the  r and  z directions.  A sphere  with  a 
radius  of  35  cells  (1.75  x 10^^  cm)  and  constant  density  1.0  x 10~^®  ^cm“^  is  placed 
on  the  grid.  Figure  (2.6)  shows  the  shape  of  the  initial  density  distribution.  The 
shape  is  not  truly  spherical  since  the  individual  cells  are  all  rectangular.  The  density 
outside  the  sphere  is  1.0  x 10^^^  gcm~^.  A constant  temperature  oi  5 K is  applied 
throughout  the  entire  grid,  and  an  ideal  gas  equation  of  state  is  used  with  7=  1.4.  For 
the  temperature  and  interior  density  chosen,  the  Jeans  length  is  less  than  the  radius. 
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Figure  2.6:  Initial  density  distribution  for  the  collapsing  sphere.  The  sphere  has 
radius  a — 1.75  x 10^^  cm  and  interior  density  p*  = 1.0  x 10~^^  gcm~^.  The  exterior 
density  is  pe  = 1.0  x 10“^°  pcm“^. 
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Aj  = 1.17  X 10^^  cm  = 23  A cells.  The  numerical  and  analytical  solutions  for  the  poten- 
tial at  t = 0 are  shown  in  figure  (2.7)  along  with  the  difference  and  relative  difference 
between  them.  The  numerical  and  analytic  solutions  show  good  agreement  consider- 
ing that  the  numerical  approximation  is  not  a true  sphere.  The  numerical  solution 
for  the  potential  is  everywhere  greater  than  the  corresponding  analytic  solution.  The 
maximum  difference  of  5.4%  occurs  at  the  edge  of  the  sphere. 

A sphere  exceeding  the  critical  Jeans  length  will  collapse  in  roughly  a free-fall 
time  given  by 

For  the  sphere  considered  here,  iff  = 2.10  x 10^^  s.  The  collapse  of  the  sphere  is 
depicted  in  figure  (2.8)  with  density  contours  displayed  at  different  fractions  of  the 
free-fall  time.  At  t = tff,  the  density  at  the  center  of  the  collapsed  sphere  has  reached 
a value  of  4.88  x 10~^^  gcm~^,  which  is  nearly  the  density  needed  for  a single  cell 
to  be  Jeans  unstable.  The  collapse  cannot  be  followed  further  without  finer  mesh 
resolution. 

2.3  Cooling  Routines 

Thermodynamic  processes  play  a key  role  in  the  evolution  of  interstellar  ma- 
terial. Indeed,  an  efficient  cooling  mechanism  is  vital  to  the  process  of  star  formation. 
In  dense,  cool  molecular  clouds,  the  dominant  form  of  cooling  is  the  collisional  ex- 
citation of  atoms  and  molecules  followed  by  the  spontaneous  emission  of  photons. 
The  dominant  atomic  species  H and  He  lack  the  low  energy  transitions  to  be  efficient 
coolants  at  low  temperatures.  The  most  abundant  molecule,  H2,  is  an  ineffective 
coolant  at  low  temperatures  due  to  a lack  of  low-lying  rotational  energy  states.  At 
high  densities,  H2  cooling  is  strongly  quenched  by  collisional  de-excitation  because 
the  molecule  lacks  an  electric  dipole  moment.  Therefore,  carbon  and  oxygen  bearing 
molecules  play  a large  role  in  driving  the  gas  cooling  in  molecular  clouds.  For  the 
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Figure  2.7:  Comparison  between  numerical  and  analytical  solutions  for  the  potential 
of  a solid  sphere,  a)  results  from  the  numerical  algorithm;  b)  results  from  the  analytic 
formulae;  c)  (numerical-analytic)  d)  (numerical-analytic)/analytic. 
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Figure  2.8:  Time  evolution  of  a collapsing  sphere.  The  plots  show  the  density  contours 
at  times  a)0.4ty^;  b)0.6t//;  c)0.8t//;  andd)1.0t//.  The  free-fall  time  t/y  = 2. lOx  10^^  s. 
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low  temperatures  and  densities  considered  in  this  study,  CO  is  the  major  coolant, 
with  O playing  a larger  role  as  the  temperature  increases  and  H2O  contributing  in 
the  shocked,  heated  regions. 

For  a gas  with  density  p and  pressure  p,  the  First  Law  of  Thermodynamics 
can  be  written; 

I -?i 

The  term  {H  — C)  accounts  for  the  change  in  heat,  (H  — C)  — pdQ/dt,  where  H is  the 
heating  rate  and  C is  the  cooling  rate.  For  a gas  obeying  the  ideal  gas  law,  equation 
(2.9),  the  First  Law  can  be  recast  into  an  equation  for  the  temperature  evolution. 


dT  . .Tdp  T.H 
~di  ~ ~ ^ T}~C 


!)• 


(2.43) 


The  cooling  time-scale  tc  is  defined  as  tc  = p/[{j  — 1)0*].  For  the  case  in  which 
the  cooling  rate  is  much  greater  than  the  heating  rate  {H/C<^1),  equation  (2.43) 
reduces  to 


dt  p dt  tc 


(2.44) 


The  first  term  on  the  right  hand  side  of  equation  (2.44)  represents  the  tempera- 
ture change  due  to  adiabatic  processes.  The  second  term  represents  the  change  in 
temperature  due  to  radiative  cooling.  The  cooling  time-scale  is  much  less  than  the 
dynamical  time-scale  of  the  problem.  The  temperature  is  updated  in  each  cell  once 
per  timestep,  and  a minimum  temperature  (typically  5 K)  is  maintained  to  account 
for  diffuse  heating  by  cosmic  rays.  The  cooling  is  applied  during  the  Lagrangian 
phase  of  the  ICE-ALE  cycle. 

The  cooling  rates  used  are  the  tabulated  rates  of  Neufeld,  Lepp,  and  Melnick 
(1995,  hereafter  NLM).  The  models  of  NLM  provide  the  total  cooling  rates  for  dense 
molecular  gas  over  a wide  range  of  density  and  temperature.  The  rates  include  a full 
consideration  of  chemical  abundances  for  389  species  made  up  from  the  elements  H, 
He,  O,  C,  N,  Si,  S,  Na,  Cl,  Fe,  P,  and  Mg.  Cooling  functions  were  obtained  for  H2, 
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Figure  2.9:  Cooling  rates.  This  plot  is  a reproduction  of  figure  (3a)  from  Neufeld  et 
al.  1995.  Each  curve  corresponds  to  a different  temperature,  with  T = \0K  being  the 
lowest  curve  and  T = 10^  '*  K the  highest. 


H2O,  CO,  O2,  HCl,  C,  and  O.  The  total  cooling  rates  are  displayed  in  figure  (2.9), 
which  is  a reproduction  of  figure  (3a)  of  NLM  extended  to  lower  densities.  Each 
curve  tracks  the  cooling  power  versus  H2  number  density  for  a constant  temperature. 
The  temperatures  range  from  T=10K  (the  lowest  curve)  to  T—IO^'^K  (the  highest 
curve)  in  logarithmic  intervals  of  0.1  dex  in  temperature.  A large  increase  occurs 
at  high  densities  for  T>300  K due  to  an  increase  in  H2O  abundance  and  the  large 
dipole  moment  of  the  water  molecule.  For  a given  density  and  temperature,  a cooling 
rate  is  found  by  linear  interpolation  between  tabulated  data  values.  In  determining 
cooling  rates,  NLM  define  an  optical  depth  parameter  in  terms  of  the  local  velocity 
gradient  dv^/dz  as 


N{M)  = 


Fn(M) 


(2.45) 


\dvz!dz\  ’ 

where  n(M)  is  the  number  density  of  species  M,  and  F is  a geometric  form  factor  set 
to  unity  for  the  plane  parallel  case.  A value  of  N{M)  = 10^°  cm“^  s is  used  in 
this  study  because  it  is  the  lowest  optical  depth  available. 
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Although  no  truly  comprehensive  model  for  interstellar  cooling  has  been  de- 
veloped, the  NLM  models  are  considered  to  be  among  the  best.  The  NLM  models 
include  the  effects  of  many  chemical  species  in  determining  the  cooling  rates,  but  the 
models  do  have  some  limitations  and  possible  sources  of  error.  The  models  consider 
only  steady  state  chemistry  and  ignore  time  dependent  processes  such  as  the  evapo- 
ration of  icy  grain  mantles.  Direct  cooling  by  grains  is  also  not  included.  Evidence 
exists  which  suggests  that  dust  cooling  dominates  behind  shocks  in  molecular  clouds 
(Whitworth  and  Clarke  1997).  If  this  is  the  case,  the  NLM  models  underestimate  the 
cooling  rates  in  shocked  regions.  NLM  assume  that  heating  due  to  an  infrared  dust 
continuum  or  radioactive  decay  is  negligible,  and  they  only  consider  fully-shielded  re- 
gions. If  external  ultraviolet  radiation  is  able  to  penetrate  a region,  photodissociation 
will  alter  the  chemical  abundances  in  the  cloud.  Recent  observations  of  the  Orion 
BN/KL  region  may  indicate  that  the  NLM  models  overestimate  the  abundances  of 
the  two  important  molecular  species  H2O  and  O2  (Bergin  et  al.  2000).  Although 
CO  is  the  dominant  coolant  for  low  temperature,  low  density  gas,  the  NLM  model 
predicts  that  O2  should  account  for  about  20%  of  the  cooling  at  low  density.  Water 
contributes  about  10%  of  the  cooling  power  at  low  temperatures.  If  the  observations 
are  typical  of  molecular  clouds,  the  NLM  models  overestimate  the  overall  cooling 
rates  by  a nontrivial  amount  for  regions  with  low  temperatures  or  low  densities.  In 
shocked  regions,  the  abundance  of  water  is  observed  to  be  much  higher  than  in  quies- 
cent regions.  So  water  remains  an  important  coolant  in  shocked  regions,  as  predicted 
by  NLM.  The  possible  sources  of  error  in  the  NLM  models  are  offsetting  to  some 
degree.  The  effects  of  altering  the  cooling  rates  are  investigated  by  parameterizing 
the  amount  of  cooling  for  a sequence  of  otherwise  identical  models.  A coefficient  is 
applied  to  the  cooling  term  in  equation  (2.44)  to  allow  various  degrees  of  cooling  from 
zero  to  the  full  NLM  cooling  rates.  The  results  are  presented  in  a later  chapter. 


CHAPTER  3 

SELF-GRAVITY  DRIVEN  INTERFACIAL  INSTABILITIES 

Interfacial  instabilities  are  an  important  feature  of  the  ISM.  Dynamic  gas  flows 
driven  by  a variety  of  mechanisms  constantly  interact  in  interesting  ways.  For  exam- 
ple, a velocity  shear  at  the  interface  between  two  fluids  drives  the  Kelvin-Helmholtz 
instability.  The  Rayleigh- Taylor  instability  arises  when  a fluid  is  accelerated  into 
a more  tenuous  fluid.  In  a recent  paper,  Hunter,  Whitaker,  and  Lovelace  (1997, 
hereafter  HWL)  identified  a new  type  of  instability  which  acts  upon  an  interface  of 
discontinuous  density  and  is  driven  by  self-gravity.  In  this  chapter,  I discuss  this 
self-gravity  interfacial  instability  (SGI)  in  Cartesian  geometry.  The  theory  of  HWL 
is  summarized,  the  computational  models  are  described,  and  results  are  presented. 

3.1  Theory 

The  discovery  of  the  SGI  came  about  as  part  of  an  ongoing  study  of  interfacial 
instabilities  in  the  ISM.  To  explore  the  nature  of  interfacial  instabilities.  Hunter  and 
Whitaker  (1987,  1989)  considered  the  stability  of  two  semi-infinite,  uniform  fluids 
in  relative  motion  having  constant,  but  generally  unequal,  unperturbed  densities, 
temperatures,  and  magnetic  held  strengths.  They  used  a normal  mode  analysis  to 
study  compressible  Kelvin-Helmholtz  instabilities  and  identify  a largely  supersonic 
thermal-dynamic  instability  which  arises  due  to  a coupling  between  the  dynamics 
and  thermodynamics  of  the  fluids.  As  a next  step,  HWL  included  self-gravity  in  the 
analysis.  By  doing  so,  they  were  able  to  identify  the  SGI,  an  interfacial  instability 
driven  by  self-gravity  that  acts  upon  any  interface  of  discontinuous  density.  The  SGI 
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persists  in  the  static  limit  and  has  a characteristic  growth  rate  of  the  order  of  the 
free-fall  time  in  the  denser  fluid. 

In  their  initial  analysis  of  the  SGI,  HWL  considered  non-stratified,  compress- 
ible media  with  an  interface  characterized  by  a velocity  slip  as  well  as  a density 
discontinuity.  In  the  analysis  below,  I consider  only  the  static,  incompressible  case. 
A full  normal  mode  analysis  for  the  determination  of  the  growth  rate  is  presented  by 
HWL  and  not  reproduced  here.  HWL  also  present  a truncated  derivation  of  the  SGI 
growth  rate  using  the  energy  principle  method  of  Bernstein  et  al.  (1958).  The  result 
is  the  same  in  the  static,  incompressible  limit.  A planar  interface  at  y = 0 is  placed  in 
equilibrium  between  two  regions  of  different  densities,  p\  below  the  interface  and  p2 
above.  The  background  acceleration  (effective  gravity)  \s  g = gy,  where  y is  a unit 
vector  in  the  y-direction.  A distinction  must  be  made  between  g and  the  accelera- 
tion due  to  the  self- gravitational  potential  <l>.  A linear  perturbation  is  applied  to  the 
equilibrium  configuration  such  that  any  perturbed  quantity  rf  is  written 

= (3.1) 


The  perturbation  results  in  the  linear  equations. 


dt 


-V  • (pu'), 


= -Vp'  + gp'  + pV$'  + [(V$)p'], 
dp' 

— = —v'  ■ Vp  — 7pV  • v'. 


(3.2) 

(3.3) 

(3.4) 


In  the  above,  p is  the  fluid  pressure,  v the  velocity,  and  7 the  ratio  of  specific  heats. 
The  perturbed  Poisson  equation  is 


= -47tGp'.  (3.5) 

In  the  incompressible  limit,  p'  = 0 except  at  the  interface,  so  the  bracketed  term  in 
equation  (3.3)  can  be  dropped.  To  examine  the  perturbed  equations,  a first  order. 
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Lagrangian  displacement  e is  defined  such  that  the  perturbed  velocity  is  given  by 
v'  — de/dt  (Bernstein  et  al.  1958).  The  equations  for  the  perturbed  density  and 
pressure  can  then  be  expressed  as 

p'  = -V  ■ (pe)  (3.6) 


and 


p = —{e  ■ V)p  — 7pV  • e. 


(3.7) 


The  perturbed  momentum  equation  takes  the  form 


(3.8) 


where  F is  a linear,  self-adjoint  operator  on  e.  For  a perturbation  of  the  form  given 
in  equation  (3.1),  the  above  equation  becomes 


(j?pe  --  — F(e). 


(3.9) 


Upon  forming  the  dot  product  of  this  equation  with  e*  (the  complex  conjugate  of  e), 
integration  over  space  and  integration  by  parts  leads  to  an  energy  principle. 


UJ 


J(fxpe*-e=  Jd^x[yp  iV-ej^  -f-(V-e:*)(e-V)p]-l-  J(fxg  e*[V-{p£)]  — J(fxp£*-V^  . 


(3.10) 

The  integrals  appearing  on  the  right  hand  side  of  equation  (3.10)  have  evident  physical 
interpretations:  the  first  describes  the  generally  stabilizing  infiuence  of  compressibil- 
ity, the  second  is  associated  with  the  familiar  Rayleigh-Taylor  (RT)  instability,  and 
the  third  represents  the  destabilizing  effect  of  self-gravity  upon  the  interface.  For  the 
incompressible  case,  V-e  = V-e*  = 0,  so  the  first  two  terms  drop  out.  The  remaining 
integrals  are  evaluated  to  give  the  incompressible  growth  rate, 

^2  _ ~2ttG{p2  - Pi)^  ^ gk{p2  - Pi) 

P2  + Pi  P2  + Pi 

The  second  term  in  equation  (3.11)  is  the  growth  rate  for  the  incompressible  RT  in- 
stability. The  first  term  gives  the  new  result  of  HWL,  the  incompressible  growth  rate 
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for  the  SGI.  Several  important  differences  are  apparent.  Self-gravity  is  destabilizing 
across  any  interface  of  discontinuous  density,  whereas  the  interface  is  RT  unstable 
only  if  g{p2  — Pi)  < 0.  For  the  SGI,  the  growth  rate  depends  upon  the  density  of 
the  environment  as  well  as  the  density  contrast  at  the  interface,  but  not  upon  the 
perturbation  wavelength.  In  contrast,  the  growth  of  the  RT  instability  depends  up)- 
on  the  density  contrast  and  the  perturbation  wavelength,  but  not  the  background 
density.  For  a fixed  value  of  g,  the  SGI  grows  more  rapidly  than  the  RT  instability 
in  relatively  dense  regions  with  relatively  long  wavelength  perturbations,  such  that 
A = 2n/k  > g/G\p2  - pi\. 

The  physical  origin  of  the  SGI  can  be  explained  as  follows  (HWL).  For  speci- 
ficity, let  the  density  below  the  interface  (pi)  be  greater  than  the  density  above  the 
interface  (p2).  The  first  term  in  equation  (3.11)  can  be  traced  to  the  pressure  d- 
ifference  across  the  interface.  An  upward  distortion  of  the  interface,  £y(x),  results 
in  a mass  excess  proportional  to  (pi  — p2)£y{x).  This  in  turn  produces  a positive 
perturbation  of  the  gravitational  potential  thus,  $'  > 0 in  the  nearby  space  above 
and  below  the  interface.  An  increase  in  the  fluid  pressure  results  owing  to  Bernoulli’s 
equation,  which  in  the  static  limit  reduces  to  p = p^' . The  potential  is  continuous 
across  the  interface,  so  the  pressure  below  the  interface,  p\  — P\^' , is  larger  than  that 
above,  p'2  = P2^  • The  pressure  difference  across  the  interface  causes  the  distortion 
to  grow. 

For  their  primary  analysis,  HWL  considered  non-stratified  media  above  and 
below  the  interface,  such  that  the  density  is  constant  in  each  region.  In  addition,  they 
considered  an  isothermal  case  with  density  stratifications  deduced  from  the  solution 
to  Poisson’s  equation.  They  found  that  the  growth  rate  for  the  stratified  case  has 
the  same  form  as  equation  (3.11)  provided  the  perturbation  wavelength  is  much 
smaller  than  the  Jeans  length.  Therefore,  numerical  models  with  stratified  density 
distributions  are  expected  to  show  SGI  growth  rates  in  good  agreement  with  equation 
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(3.11)  when  perturbation  wavelengths  are  small,  but  the  deviation  from  equation 
(3.11)  will  grow  larger  as  the  perturbation  wavelength  increases  to  an  appreciable 
fraction  of  the  Jeans  length. 

3.2  Models 

In  order  to  study  the  SGI,  models  are  given  initial  conditions  that  correspond 
to  hydrostatic  equilibrium.  Because  the  ISM  is  a turbulent  medium,  a configuration  of 
hydrostatic  equilibrium  is  extremely  unlikely  to  occur,  so  the  results  of  these  models 
are  not  directly  comparable  to  observations.  However,  an  equilibrium  is  necessary  in 
order  to  study  the  SGI  in  isolation  without  other  processes  interfering.  Once  this  is 
accomplished,  more  comprehensive  models  of  the  ISM  can  be  formed  to  study  the 
effects  of  the  SGI  when  coupled  to  other  phenomena,  such  as  magnetic,  thermal,  or 
dynamic  instabilities.  For  instance,  velocity  shear  must  be  added  to  test  the  theory  of 
HWL  that  the  SGI  enhances  the  growth  of  the  Kelvin-Helmholtz  instability  for  small 
to  moderate  Mach  numbers.  This  type  of  problem  is  left  for  later  work.  All  models  in 
this  chapter  are  presented  with  time  units  in  terms  of  the  theoretical  e-folding  time 
te,  defined  as  the  time  needed  for  a perturbed  quantity  to  grow  by  a factor  of  e.  The 
e-folding  time  is  the  reciprocal  of  the  theoretical  growth  rate,  Itg  = 

When  a fluid  is  in  hydrostatic  equilibrium,  the  forces  due  to  self-gravity  and 
thermal  pressure  exactly  balance,  so  the  gas  remains  at  rest.  Typically,  hydrostatic 
equilibrium  is  invoked  in  spherical  models  of  stars.  However,  the  same  principle  may 
be  applied  in  other  geometries  as  well.  The  SGI  models  employ  a square  Gartesian 
{x,  y)  grid  with  either  129  or  257  cells  in  each  direction.  The  2-D  equilibrium  problem 
with  planar  symmetry  reduces  to  a 1-D  problem;  only  the  direction  normal  to  the 
interface  need  be  considered.  The  equilibrium  configuration  for  each  region  is  solved 
for  independently  from  the  other  as  a function  of  y,  where  y>0  is  the  distance  from 
the  interface.  The  initial  density  and  temperature  distributions  are  determined  by  si- 
multaneously solving  the  equation  of  hydrostatic  equilibrium  and  Poisson’s  equation. 
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The  equation  of  hydrostatic  equilibrium  is 

1 dp 

pdy  dy' 

A polytropic  equation  of  state  is  used,  so  that, 

p = Kp'^, 


(3.12) 


(3.13) 


where  the  constant  k is  determined  by  the  values  of  p and  p at  the  interface,  n = 
Po/Po^-  The  poly  tropic  relation  allows  the  replacement  of  the  pressure  in  the  hydro- 
static equation  with  a simple  function  of  the  density.  Poisson’s  equation,  equation 
(2.34),  simplifies  in  the  1-D  Cartesian  case  to 


dy"^ 


= — 47rGp. 


(3.14) 


The  combination  of  equations  (3.12-3.14)  results  in  a differential  equation  for  the 
density  distribution. 


d^P  ^ (7-2)  ^ 47rGp(^ 

dy'^  p \dy)  7« 


= 0. 


(3.15) 


For  the  special  case  that  7 = 2,  the  second  term  drops  out  leaving 

(fp  2nGp 
dy"^ 


K 


(3.16) 


The  boundary  values  for  this  problem  are  such  that  at  the  interface,  the  density  is  a 
maximum  (p  = po)  and  the  gradient  of  the  density  goes  to  zero  {dp/dy  =0).  Under 
these  conditions,  the  solution  to  equation  (3.16)  is 


p = pocos^^).  (3.17) 

Valid  densities  are  found  from  ?/  = 0 to  {y/a)  = tt/2,  although  for  best  results  the 
scale  height  a is  chosen  so  that  y/a  is  still  significantly  less  than  tt/2  at  the  grid 
boundaries.  The  scale  height  is  determined  by  the  density  and  temperature  at  the 
interface.  For  7 = 2, 

2 K \ RTq 
2nG  2ttG  p po  ’ 


(3.18) 
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where  R and  //  are  the  gas  constant  and  molecular  weight,  respectively.  The  ideal  gas 
equation  of  state,  equation  (2.9),  can  be  used  to  determine  the  temperature  distri- 
bution once  the  density  and  pressure  are  known.  For  the  7 = 2 case,  the  temperature 
follows  the  same  distribution  as  the  density, 

T = ToCos^^).  (3.19) 

Although  unphysical  and  not  applicable  to  the  ISM,  the  choice  of  7 = 2 is  useful 
for  several  reasons.  The  very  stiff  equation  of  state  aids  in  limiting  models  to  the 
incompressible  regime.  More  importantly,  the  analytic  solution  was  used  in  debugging 
the  computational  method  used  for  all  SGI  models.  The  analytic  solution  provides  a 
means  by  which  to  test  the  numerical  solver  developed  for  general  values  of  7. 

In  order  to  better  represent  the  ISM,  a method  for  determining  equilibrium 
conditions  for  adiabatic  models  with  any  value  of  7 has  been  developed.  The  general 
solution  of  equation  (3.15)  requires  a numerical  integration.  The  problem  is  simplified 
by  recasting  the  problem  in  the  so-called  Emden  variables  (e.g.  Chandrasekhar  1939). 
A dimensionless  variable  9 is  used  to  parameterize  the  temperature,  density,  and 
pressure  of  a polytropic  gas  as 

T = Tq9  p = Po9^  p = (3.20) 


As  before,  the  subscripts  denote  values  at  the  interface.  The  polytropic  index  n is 
related  to  7 by  n = 1/(7  — 1).  Rewriting  equation  (3.15)  as  a function  of  9 results  in 
a simpler  equation, 

= -9^.  (3.21) 

dy^ 

As  in  the  7 = 2 case,  the  scale  height  a is  determined  from  the  density  and  temperature 
at  the  interface.  For  general  n. 


a 


2 


«(n  -I-  1) 


Po 


( 


1 — n 
n 


) 


(n  1 ) i?  Tq 
47tG  p Po 


AttG 


(3.22) 
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The  problem  is  simplified  further  by  the  definition  of  another  dimensionless  variable  ^ 
such  that  y = a^.  When  ^ is  substituted  into  equation  (3.21),  the  result  is  a Cartesian 
analog  to  the  Lane-Emden  equation, 


= -r. 


(3.23) 


The  choice  of  7 = 2 (n  = 1)  reproduces  equation  (3.16).  A particularly  simple  solution 
to  equation  (3.23)  exists  for  the  strictly  incompressible  case  7 = oo  (n  = 0).  As 
mentioned  in  HWL,  the  isothermal  case  7 = 1 (n  = 00)  has  an  analytic  solution,  as 
does  the  7 = 2/3  (n=  — 3)  case.  To  solve  for  the  equilibrium  state  for  any  other  value 
of  7,  a fourth  order  Runge-Kutta  numerical  integration  scheme  (Press  et  al.  1992)  is 
employed.  The  numerical  and  analytic  solutions  show  excellent  agreement  for  both 
the  7 = 1 and  7 = 2 cases.  As  an  example.  Table  (3.1)  compares  the  numerical  values 
of  ^(^)  with  the  analytical  values  of  cos(^)  in  the  denser  region  for  the  7 = 2 model 
masgl29dl.  The  agreement  is  perfect  to  the  edge  of  the  grid.  The  Runge-Kutta 
solver  uses  a resolution  four  times  finer  than  the  CFDLIB  grid  for  improved  accuracy 
and  to  obtain  0(^)  and  d9/d^  at  both  cell  centers  and  cell  faces.  Having  determined 
0(^),  density  and  temperature  values  throughout  the  grid  are  assigned  using  equations 
(3.20).  An  example  of  the  initial  density  and  temperature  distributions  for  the  SGI 
model  cnsg257el  is  shown  in  figure  (3.1).  This  model  has  a relatively  small  scale 
height,  which  results  in  the  steep  density  decline  in  the  denser  medium. 

In  setting  up  an  equilibrium  configuration,  boundary  conditions  must  be  care- 
fully considered.  For  constraining  the  dynamic  processes  of  the  simulation  within 
CFDLIB,  reflective  boundary  conditions  have  proven  most  successful  in  maintaining 
pressure  equilibrium  along  all  four  edges.  Two  types  of  boundary  conditions  are  im- 
plemented in  MUDPACK  to  confine  the  potential  solution.  Along  the  left  and  right 
boundaries,  periodic  boundary  conditions  are  used.  This  is  a logical  choice  given  the 
symmetry  of  the  problem  and  results  in  the  necessity  for  the  perturbation  function 
to  be  periodic  as  well.  The  top  and  bottom  boundaries  are  constrained  by  values  of 
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Table  3.1:  Comparison  of  numerical  and  analytical  results 


i = y/a 

numerical  0{^) 

analytical  cos(^) 

l.E+16 

0.0039269908 

0.999992289 

0.999992289 

5.E+16 

0.019634954 

0.99980724 

0.99980724 

9.E+16 

0.0353429172 

0.999375504 

0.999375504 

1.3E+17 

0.0510508804 

0.998697187 

0.998697187 

1.7E+17 

0.0667588436 

0.997772456 

0.997772456 

2.1E+17 

0.0824668068 

0.99660154 

0.99660154 

2.5E+17 

0.09817477 

0.995184727 

0.995184727 

2.9E+17 

0.113882733 

0.993522367 

0.993522367 

3.3E+17 

0.129590696 

0.99161487 

0.99161487 

3.7E+17 

0.14529866 

0.989462708 

0.989462708 

4.1E+17 

0.161006623 

0.98706641 

0.98706641 

4.5E+17 

0.176714586 

0.984426568 

0.984426568 

4.9E+17 

0.192422549 

0.981543834 

0.981543834 

5.3E+17 

0.208130512 

0.978418919 

0.978418919 

5.7E+17 

0.223838476 

0.975052593 

0.975052593 

6.1E+17 

0.239546439 

0.971445688 

0.971445688 

6.5E+17 

0.255254402 

0.967599093 

0.967599093 

6.9E+17 

0.270962365 

0.963513757 

0.963513757 

7.3E+17 

0.286670328 

0.959190689 

0.959190689 

7.7E+17 

0.302378292 

0.954630955 

0.954630955 

8.1E+17 

0.318086255 

0.94983568 

0.94983568 

8.5E+17 

0.333794218 

0.944806047 

0.944806047 

8.9E+17 

0.349502181 

0.939543297 

0.939543297 

9.3E+17 

0.365210144 

0.93404873 

0.93404873 

9.7E+17 

0.380918108 

0.928323699 

0.928323699 

l.OlE+18 

0.396626071 

0.922369619 

0.922369619 

1.05E+18 

0.412334034 

0.916187958 

0.916187958 

1.09E+18 

0.428041997 

0.909780241 

0.909780241 

1.13E+18 

0.44374996 

0.903148049 

0.903148049 

1.17E+18 

0.459457924 

0.89629302 

0.89629302 

1.21E+18 

0.475165887 

0.889216843 

0.889216843 

1.25E+18 

0.49087385 

0.881921265 

0.881921265 

1.29E+18 

0.506581813 

0.874408087 

0.874408087 

Tn  this  notation,  5.E+16  = 5.0  x 10^®  cm. 
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Y cell 


Figure  3.1;  Initial  density  and  temperature  distributions  for  SGI  model  cnsg257el. 
The  temperature  axis  has  been  displaced  to  better  distinguish  between  the  density 
and  temperature. 


the  potential  gradient  d^/dy  determined  for  the  initial  equilibrium  distribution.  The 
physical  derivative  is  obtained  from  the  dimensionless  derivative  via 


^ _ (ra  + \)RTq  fM 
dy  y,a  \d^ 


(3.24) 


The  gradients  are  calculated  before  the  simulation  begins.  In  perturbed  models, 
the  validity  of  the  gradient  values  at  the  top  and  bottom  boundaries  decreases  as  the 
perturbations  grow.  However,  since  the  total  mass  within  the  grid  is  held  constant  by 
the  reflective  boundary  conditions,  the  spurious  motions  induced  by  not  considering 
the  time  dependence  of  the  gradient  values  at  the  boundaries  are  very  small  compared 
to  the  motions  of  interest  near  the  interface.  The  gravity  solver  begins  to  fail  only 
when  a perturbation  grows  large  enough  to  interact  directly  with  the  boundaries. 

The  discussion  above  explains  the  methods  used  to  create  a state  of  hydrostatic 
equilibrium  for  a self-gravitating  gas.  As  I show  later,  the  equilibria  achieved  are  very 
good.  However,  representing  a discontinuity  on  a discrete  grid  with  finite  cell  sizes 
cannot  be  done  without  introducing  some  numerical  artifacts  into  the  solution.  A 
perfect  equilibrium  is  not  attainable;  one  can  only  be  approximated.  Setting  up 
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the  problem  on  a computational  grid  requires  careful  consideration  of  whether  to 
place  the  interface  at  a cell  center  or  a cell  face.  CFDLIB  uses  cell-centered  state 
variables,  so  the  density,  temperature,  pressure,  and  velocity  terms  are  initialized  at 
cell  centers  to  begin  the  simulation.  However,  face-centered  velocity  and  pressure 
terms  are  also  used  in  the  computational  cycle.  Experiments  indicate  that  a better 
equilibrium  is  achieved  when  the  interface  is  placed  at  a cell  face.  At  an  interface 
of  discontinuous  density,  pressure  balance  between  two  perfect  gases  with  the  same 
7 and  molecular  weight  is  achieved  by  altering  the  temperatures  so  the  product  of 
temperature  and  density  remains  constant.  The  different  density  and  temperature 
values  on  either  side  of  the  grid  result  in  different  scale  heights,  as  described  by 
equation  (3.22).  As  a consequence  of  placing  the  interface,  and  hence  po  and  Tq,  at 
a cell  face,  the  maximum  cell-centered  values  in  either  region  are  removed  from  po 
and  To  by  half  a cell.  Likewise,  the  maximum  cell-centered  pressure  is  slightly  lower 
than  po  in  both  regions.  Since  a configuration  of  discontinuous  density  results  in  a 
different  scale  height  on  either  side  of  the  interface,  the  pressure  falls  off  at  different 
rates  in  the  two  regions.  As  a result,  the  simulation  begins  with  a slight  pressure 
imbalance  at  the  interface.  The  more  tenuous  fluid  has  a larger  scale  height,  so  the 
cell-centered  pressure  on  that  side  of  the  interface  is  higher  than  the  pressure  on  the 
side  containing  the  denser  fluid.  If  the  scale  heights  are  increased  in  both  regions,  the 
pressure  difference  across  the  interface  decreases.  A larger  scale  height  also  results  in 
smaller  gradients  for  the  pressure,  density,  temperature,  and  potential,  which  helps 
to  maintain  a stable  equilibrium.  If  the  only  goal  were  to  create  the  best  equilibrium 
possible,  a very  large  scale  height  would  be  desired.  However,  a very  large  scale 
height  is  not  always  obtainable  when  performing  simulations  on  astrophysical  scales. 
The  effects  of  having  relatively  short  scale  heights  are  examined  on  static  equilibria 
before  considering  perturbed  cases. 
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Although  an  in-depth  consideration  of  the  RT  instability  is  not  the  goal  of 
this  study,  it  is  of  interest  to  compare  the  morphological  similarities  and  differences 
between  the  SGI  and  the  RT  instability.  To  provide  a better  basis  for  comparisons, 
RT  models  are  also  begun  from  an  initial  state  of  hydrostatic  equilibrium.  This 
requires  a more  careful  construction  than  that  in  the  example  of  chapter  2.  The 
underlying  method  is  the  same  as  that  used  for  self-gravitating  models,  but  the 
equations  are  easier  to  solve.  For  a gas  that  is  not  self-gravitating,  but  feels  a constant 
acceleration  g (assumed  to  be  positive  in  the  upward  direction),  the  equation  of 
hydrostatic  equilibrium  is  simply 


dp 

T =9P- 
dy 


(3.25) 


The  gravitational  acceleration  g is  determined  for  each  model  by  setting  the  RT 
growth  rate  (the  second  term  on  the  right  hand  side  of  equation  (3.11))  equal  to  the 
SGI  growth  rate  (the  first  term  on  the  right  hand  side  of  equation  (3.11))  so  that  for 
a perturbation  wavelength  A, 


g = G{p,  - p2)X.  (3.26) 

As  a result,  SGI  and  RT  models  with  the  same  physical  parameters  at  the  interface 
share  the  same  theoretical  incompressible  growth  rate  and  e-folding  time.  Again 
assuming  a polytropic  equation  of  state,  the  equation  of  hydrostatic  equilibrium  be- 
comes 

pG-^)dp  - ^dy.  (3.27) 

7K 

Integrating  and  setting  p = po  at  the  interface  gives 
The  temperature  distribution  is 

T=^^r:^y  + To 


(3.29) 
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Figure  3.2:  Initial  density  and  temperature  distributions  for  RT  model  cart257el. 
The  temperature  axis  has  been  displaced  to  better  distinguish  between  the  density 
and  temperature. 


In  the  absence  of  self-gravity,  the  analytic  solution  applies  to  all  values  of  7,  so  no 
Runge-Kutta  solver  is  needed.  The  static  equilibria  are  typically  more  robust  for 

the  RT  models  than  for  the  SGI  models.  The  MUDPACK  routines  are  not  used 

« 

when  running  RT  models,  so  only  the  reflective  boundary  conditions  of  CFDLIB 
are  applied.  Figure  (3.2)  shows  an  example  of  the  initial  density  and  temperature 
distributions  for  the  RT  model  cart257el.  Other  than  at  the  interface,  both  density 
and  temperature  increase  in  the  direction  of  y. 

3.3  Results  for  Unperturbed  Models 

For  every  perturbed  model  considered,  a non-perturbed  model  was  run  to  test 
the  robustness  of  the  initial  equilibrium.  The  degree  to  which  equilibrium  is  main- 
tained varies  from  model  to  model,  with  generally  good  results.  In  all  descriptions 
that  follow,  I refer  to  the  denser  region  below  the  interface  and  the  more  tenuous 
region  above  the  interface  as  region  1 and  region  2,  respectively,  with  subscripts  de- 
noting values  at  the  interface.  A key  factor  is  the  relative  size  of  the  scale  height 
versus  the  grid  size.  Equilibrium  is  better  maintained  when  the  scale  height  of  the 
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Figure  3.3:  Pressure  comparison  for  model  cnsg257d0  at  f = 0 and  t = lOte-  Values 
are  plotted  along  a constant  value  of  x for  a)the  pressure  distribution  at  t = 0;  b)the 
pressure  distribution  at  t = lOtg;  c)the  difference  p(lOte)  — p(0);  and  d)the  relative 
difference  (p(lOte)  — p(0))/p(0). 


denser  fluid  ai  is  large  compared  to  the  spatial  extent  of  the  denser  region.  The 
scale  height  is  always  larger  in  region  2,  so  equilibrium  is  better  maintained  above 
the  interface  than  below.  The  quantity  Si  is  defined  as  the  ratio  of  the  scale  height 
to  the  spatial  extent  in  region  1, 


= -jr-f ^ i JT—-  (3-30) 

# oj  rows  in  region!  x cell  size 

For  example,  model  cnsg257d0  has  a scale  height  much  larger  than  the  spatial  extent 
of  region  1,  resulting  in  Sj  =3.2.  After  a long  time  (t  = lOtg),  the  maximum  y- velocity 
is  0.65 cm compared  to  a sound  speed  of  24, 124 cm This  is  effectively  a 
velocity  of  zero  since  the  gas  cannot  traverse  a single  cell  by  the  end  of  a simulation. 
The  pressure  change  over  that  time  is  very  small.  Figure  (3.3)  shows  the  pressure 
along  a line  of  constant  x for  model  cnsg257d0  at  t = 0 and  t = 10te-  The  relative  pres- 
sure difference  is  less  than  0.02%.  Model  cnsg257e0  is  identical  to  model  cnsg257d0 
in  every  way  except  that  the  cell  size  is  five  times  larger.  Because  of  the  larger  cell 
size,  the  relative  size  of  the  scale  height  compared  to  the  grid  size  is  smaller,  Si  =0.65. 
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Figure  3.4:  Pressure  comparison  for  model  cnsg257e0  at  ^ = 0 and  t = 10te. 


As  a result,  the  departure  from  equilibrium  is  higher.  The  maximum  y- velocity  at 
t = lOte  is  25.1  cm and  the  relative  pressure  change,  shown  in  figure  (3.4),  is 
about  20%  at  the  lower  boundary.  Although  the  velocity  after  lOtg  is  nearly  forty 
times  larger  in  model  cnsg257e0  than  in  model  cnsg257d0,  it  it  still  much  smaller  than 
the  sound  speed.  This  indicates  that  even  when  the  scale  height  is  small  enough  to 
result  in  very  low  density  and  temperature  values  at  the  boundaries,  the  equilibrium 
is  still  well  maintained.  Any  departure  from  equilibrium  will  not  interfere  with  the 
growth  of  induced  perturbations  at  the  interface.  The  requirement  that  the  density 
and  temperature  remain  positive  places  an  upper  limit  on  the  cell  size. 

RT  models  do  not  require  a numerical  solution  for  the  initial  equilibrium  or 
the  application  of  addition  boundary  conditions  for  the  gravitational  potential.  As  a 
result,  unperturbed  RT  models  maintain  better  static  equilibria.  Figure  (3.5)  shows 
the  pressure  along  a line  of  constant  x for  model  cart257d0  at  t = 0 and  t = lOtg-  The 
pressure  difference  is  several  orders  of  magnitude  lower  than  that  shown  in  figure  (3.3) 
for  model  cnsg257d0.  Likewise,  the  pressure  difference  for  model  cart257e0,  shown 
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Figure  3.5:  Pressure  comparison  for  model  cart257d0  at  t = 0 and  t=10te- 


in  figure  (3.6),  is  several  orders  of  magnitude  lower  than  the  pressure  difference  for 
model  cnsg257e0,  shown  in  figure  (3.4). 

3.4  Results  for  Perturbed  Models 

Perturbed  models  are  initiated  with  a velocity  perturbation  of  the  form 

v{x,y)  — Vocos{kxx)e^~’^’'^\  (3.31) 

where  Vq  is  the  maximum  velocity  at  the  interface.  The  original  set  of  models  (the 
moderate  temperature  models)  are  an  exception  in  that  a sine  function  was  used 
rather  than  a cosine  function.  The  switch  to  cosine  minimized  distortions  due  to  the 
boundary  conditions.  The  value  of  vo  is  kept  small  compared  to  the  sound  speed 
in  the  denser  medium  {vq/c\  1)  to  approximate  incompressible  behavior.  Unless 
otherwise  noted,  ky  = kx  = k.  A list  of  parameters  for  a representative  set  of  models 
is  given  in  Table  (3.2).  Although  additional  models  were  run,  specific  remarks  are 
limited  to  those  listed  in  the  table.  The  perturbation  wavelength  X = 2n/k  always 
chosen  to  be  less  than  the  Jeans  length  Xj  (see  equation  2.38).  The  ratio  J = X/Xj 
is  an  important  gauge  of  the  overall  compressibility  of  the  problem.  As  J increases. 
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Figure  3.6:  Pressure  comparison  for  model  cart257e0  at  ^ = 0 and  t = 10te. 


the  compressibility  increases.  Because  the  scale  height  and  the  Jeans  length  vary  the 
same  way  with  density  and  temperature,  a model  with  a large  J value  has  a small 
s value.  In  Table  (3.2)  only  J is  listed  due  to  its  more  significant  physical  meaning. 
The  theoretical  growth  rate  (equation  3.11)  was  developed  in  the  incompressible, 
non-stratified  limit,  which  corresponds  to  a very  low  value  for  J.  So  the  models  with 
very  low  J values  come  closest  to  matching  the  theoretical  conditions. 

A few  words  are  in  order  about  the  convention  used  in  naming  the  models. 
The  first  letter  of  each  name  describes  the  temperature  eis  moderate,  hot,  or  cold. 
The  second  letter  tells  whether  the  initial  equilibrium  is  determined  by  analytic  or 
numerical  means.  The  next  two  letters  label  the  model  as  either  sg  for  SGI  or  rt 
for  RT.  Each  model  is  run  on  a square  grid  with  either  129  or  257  cells  in  each 
direction.  Models  designated  with  a d are  low  J models,  while  an  e designates  a 
high  J model.  Finally,  a number  is  used  to  otherwise  differentiate  among  models. 
Unperturbed  models  are  designated  0,  with  perturbed  models  numbered  sequentially 
higher.  Dividing  the  models  according  to  temperature  is  a useful  way  to  discuss 
the  SGI  under  very  different  conditions.  The  first  set  of  models  are  referred  to  as 
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Table  3.2:  Model  Parameters 


model 

7 

Pi 

gem  ^ 

Ti 

K 

P2 

gcm~^ 

T2 

K 

vo/ci 

kl 

cm~^ 

ky 

cm~^ 

Ji 

Ax* 

cm 

it 

s 

masgl29d0 

2.0 

2.0 

131 

1.0 

262 

0.0 

0.0 

0.0 

— 

2.0 

8.458 

masgl29dl 

2.0 

2.0 

131 

1.0 

262 

0.01 

9.74 

9.74 

.0403 

2.0 

8.458 
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moderate  temperature  models.  Many  of  these  models  were  run  in  the  testing  phase  of 
the  code,  but  little  analysis  was  performed.  The  bulk  of  the  science  was  accomplished 
with  high  and  low  temperature  models.  Because  the  scale  height  and  Jeans  length 
scale  as  hot  models  generally  have  high  s and  low  J values.  Cold  models  have 
smaller  s and  larger  J values.  The  hot  models  come  closer  to  approximating  the 
incompressible,  non-stratified  conditions  under  which  the  theoretical  growth  rates 
hold.  The  cold  models  are  more  representative  of  the  conditions  in  cold  molecular 
clouds. 

3.4.1  Moderate  Temperature  Models 

The  first  set  of  models  is  of  moderate  temperature  (Ti  = 131  K)  and  uses  the 
analytic  7 = 2 solution  for  the  initial  equilibrium.  All  moderate  temperature  models 
have  low  J values  and  a density  ratio  of  2.  They  primarily  were  used  for  testing  the 
procedure  for  establishing  the  initial  hydrostatic  equilibrium.  An  example  is  model 
masgl29dl,  the  evolution  of  which  is  shown  in  figure  (3.7).  As  previously  noted, 
model  masgl29dl  was  perturbed  using  a sine  function  with  amplitude  .Olci  (ci  = 
104, 294  cm  s“^).  Asymmetries  arise  early  in  the  simulation  due  to  boundary  effects. 
Because  of  the  asymmetries,  little  analysis  was  done  on  the  moderate  temperature 
models.  Even  so,  model  masgl29dl  clearly  shows  one  large  difference  between  the 
SGI  and  the  RT  instability.  Whereas  the  RT  instability  is  characterized  by  dense 
columns  of  material  moving  into  and  displacing  the  tenuous  material,  the  SGI  does 
the  opposite.  The  SGI  gives  rise  to  tenuous  columns  moving  into  the  denser  material. 
This  trait  is  common  to  all  SGI  models. 

According  to  HWL,  the  SGI  growth  rate  for  incompressible  fluids  does  not 
depend  upon  the  perturbation  wavelength.  To  test  the  validity  of  this  claim,  a series 
of  simulations  was  performed  in  which  only  the  value  of  kx  was  altered.  A sample 
of  the  results  is  shown  in  figure  (3.8).  The  value  of  ky  was  kept  constant  amongst 
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Figure  3.7:  Evolution  of  the  SGI  for  model  masgl29dl.  Filled  contours  show  the 
density  distribution  at  a)lte;  b)2te;  c)3te;  and  d)4te.  The  e-folding  time  is  te  = 
8.458  X 10^^  s.  The  cell  size  Aa;  = Ay  = 2.0  x 10^®cm. 
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Figure  3.8:  Comparison  of  moderate  temperature  models  with  varying  The  pan- 
els show  the  density  contours  at  t = 3tg  for  models  a)masgl29d2;  b)meisgl29dl; 
c)mcisgl29d3;  and  d)masgl29d4.  The  SGI  grows  at  the  same  rate  in  all  four  simu- 
lations. All  four  models  share  the  same  e-folding  time,  = 8.458  x 10^^  s,  and  cell 
size,  Aa;  = Ay  = 2.0  x 10^®cm. 
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the  models  to  keep  the  velocity  perturbations  as  similar  as  possible.  The  prediction 
holds  up  very  well. 

3.4.2  High  Temperature  Models 

High  temperature  models  provide  the  best  comparison  between  the  computa- 
tional and  theoretical  results  due  to  their  low  J values.  These  models  are  character- 
ized by  temperature  Ti  = 1000  K and  density  pi  = 2 x 10“^^^  cm~^.  The  corresponding 
values  of  T2  and  p2  depend  upon  the  density  contrast  at  the  interface.  Two  cases  are 
considered  in  detail,  pi/p2  = 2 and  pi/p2  = 5.  The  sound  speed  in  the  denser  medium 
for  models  with  7 = 2 is  cj  = 288,  340  cm 

For  the  p\/ P2  — 2 case,  comparison  is  made  to  the  moderate  temperature 
model  masgl29dl.  Figure  (3.9)  shows  the  evolution  of  model  hasgl29dl.  The  only 
differences  from  model  masgl29dl  are  the  temperature,  the  perturbation  amplitude, 
and  the  use  of  a cosine  in  the  perturbation  function.  Once  again,  the  analytic  7 = 2 
solution  is  used  in  the  initialization.  The  new  perturbation  function  does  far  better 
at  preserving  the  symmetry  across  the  grid.  Otherwise,  the  SGI  morphology  and 
growth  rate  are  very  similar.  The  symmetry  of  model  hasgl29dl  helps  to  see  that 
the  initial  perturbation  wavelength  is  preserved  well  into  the  nonlinear  regime,  a trait 
the  SGI  shares  with  the  RT  instability. 

When  running  any  sort  of  computer  simulation,  it  is  necessary  to  explore  the 
effects  of  the  finite  grid  resolution.  Model  hasgl29dl  was  duplicated  on  a larger  grid 
to  produce  model  hasg257dl.  The  result  is  displayed  in  figure  (3.10).  The  cell  size 
for  model  hasg257dl  is  one-half  the  cell  size  for  model  hasgl29dl,  so  the  physical 
spatial  extent  is  about  the  same.  The  SGI  grows  at  very  nearly  the  same  rate  in 
both.  However,  the  higher  resolution  model  clearly  shows  fine  structures  not  seen  at 
lower  resolution.  This  leads  to  the  conclusion  that  although  low  resolution  is  sufficient 
for  determining  the  growth  rate  to  a good  level  of  accuracy,  a higher  resolution  is 


56 


Figure  3.9:  Evolution  of  the  SGI  for  model  hasgl29dl.  Filled  contours  show  the 
density  distribution  at  a)lte;  b)2te;  c)3tg;  and  d)4te-  The  e-folding  time  is  = 
8.458  X 10^^  s.  The  cell  size  Ax  = Ay  = 2.0  x lO^^cm. 
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Figure  3.10:  Evolution  of  the  SGI  for  model  hasg257dl.  This  model  is  a high  reso- 
lution version  of  model  hasgl29dl.  Filled  contours  show  the  density  distribution  at 
a)lte;  b)2te;  c)3te]  and  d)4te-  The  e-folding  time  is  tg  = 8.458  x 10^^  s.  The  cell  size 
Aa;  = Ay  = 1.0  x 10^®cm. 
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needed  to  examine  the  nonlinear  morphology.  I restrict  the  bulk  of  my  analysis  to 
the  consideration  of  high  resolution  models. 

In  order  to  measure  the  experimental  growth  rate  of  the  SGI,  I trace  the  ve- 
locity growth  inside  one  of  the  tenuous  columns.  For  each  model,  I plot  ln[v{t)  j v{tQ)\ 
versus  the  time  in  units  of  the  theoretical  e-folding  time.  The  slope  of  the  best-fit  line 
corresponds  to  the  growth  rate.  The  choice  of  to  is  arbitrary  since  any  change  alters 
the  line’s  intercept  but  not  its  slope.  I typically  use  to  = l^e-  Of  more  consequence 
is  the  choice  of  which  data  points  to  include  in  the  least  squares  calculation  that 
determines  the  best-fit  line.  For  some  models  the  data  are  such  that  the  inclusion  or 
exclusion  of  a few  points  can  change  the  slope  of  the  line  by  as  much  as  ten  percent. 
When  possible,  the  same  time  range  is  used  for  determining  the  growth  rates  for 
models  that  are  to  be  compared  to  one  another.  Figure  (3.11)  shows  the  result  for 
model  hasg257dl.  Multiple  modes  are  excited  by  the  perturbation,  initially  causing 
the  velocity  to  fall  due  to  the  presence  of  small  amplitude  sound  waves.  Once  the 
excited  SGI  mode  begins  to  dominate  the  motion,  a linear  relationship  is  clearly  seen. 
An  early  stage  of  exponential  growth  is  another  trait  shared  between  the  SGI  and 
the  RT  instability.  An  experimental  result  in  perfect  agreement  with  theory  would 
produce  a slope  of  unity.  For  model  hasg257dl  the  growth  rate  is  a;  = 1.028t“^,  very 
close  to  the  theoretical  prediction. 

Up  to  this  point,  all  models  have  used  the  7 = 2 analytic  formulae  to  create 
the  initial  static  equilibrium.  From  here  on,  models  with  7 = 1.4  are  considered.  The 
numerical  (Runge-Kutta)  solver  is  used  to  establish  the  initial  equilibrium.  Model 
hnsg257dl  is  a 7 = 1.4  model  which  corresponds  to  model  hasg257dl.  The  lower  7 
results  in  a lower  sound  speed  (ci  = 241, 242 cm s“^),  which  translates  into  a lower 
perturbation  amplitude  and  a higher  J value.  Otherwise,  models  hnsg257dl  and 
hasg257dl  are  identical.  As  shown  in  figure  (3.12),  the  SGI  development  is  the  same 
as  that  for  model  hasg257dl.  Figure  (3.13)  shows  model  hnsg257dl  at  a later  time 
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Figure  3.11:  Velocity  growth  for  model  hasg257dl.  Data  points  are  normalized  to 
the  velocity  value  at  to  = l.Otg-  The  best-fit  line  using  data  points  from  O.SOtg  to 
l.SOtg  yields  a growth  rate  of  1.028^7^  The  e-folding  time  is  tg  = 8.458  x 10^^  s. 


60 


Figure  3.12:  Evolution  of  the  SGI  for  model  hnsg257dl.  Filled  contours  show  the 
density  distribution  at  a)lte;  b)2tg;  c)3te;  and  d)4tg.  The  e-folding  time  is  te  — 
8.458  X 10^^  s.  The  cell  size  Ax  = Ay  = 1.0  x 10^®cm. 
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Figure  3.13:  Model  hnsg257dl  at  t = 5ig-  The  SGI  has  evolved  to  form  islands  of 
tenuous  gas  embedded  in  the  denser  gas. 
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of  t = 5te.  The  structures  seen  near  the  bottom  of  the  grid  are  due  to  the  reflective 
boundary  condition,  and  hence  are  unphysical.  Closer  to  the  interface,  islands  of 
tenuous  gas  are  seen  embedded  in  the  denser  material.  The  interface  has  moved  very 
little  into  the  tenuous  fluid,  but  the  wide  dense  bubbles  displace  a lot  of  tenuous 
fluid.  As  in  the  7 = 2 case,  the  experimental  growth  rate  of  a tenuous  column  in 
model  hnsg257dl  agrees  with  the  theoretical  growth  rate.  Figure  (3.14)  shows  the 
velocity  growth  for  model  hnsg257dl.  The  calculated  slope  gives  a growth  rate  of 
u}=  1.031^7^.  Switching  from  7==  2 to  7 = 1.4  has  a minimal  effect  on  the  growth  of 
the  SGI. 

Since  the  SGI  and  the  RT  instability  are  both  density  interface  instabilities, 
it  is  instructive  to  compare  and  contrast  the  two.  Model  hart257dl  is  an  RT  model 
with  the  same  interfacial  densities  and  temperatures  and  the  same  grid  parameters 
as  model  hnsg257dl.  From  equation  (3.26),  the  value  of  g which  gives  the  same 
theoretical  growth  rate  as  that  for  model  hnsg257dl  is  g = A.288  x 10“^^cms“^.  The 
morphology  of  model  hart257dl  is  shown  in  figure  (3.15).  In  the  early  stage  of  de- 
velopment, the  SGI  and  the  RT  instability  appear  quite  similar.  Both  preserve  the 
perturbation  wavelength  well  into  the  nonlinear  regime.  However,  at  late  times  the 
SGI  and  the  RT  instability  appear  very  different.  As  expected  for  the  RT  instability, 
widely  spaced  columns  of  dense  gas  stream  into  the  more  tenuous  medium.  (The 
characteristic  mushroom  cap  structure  arises  due  to  a secondary  Kelvin-Helmholtz 
instability.)  In  contrast,  the  SGI  gives  rise  to  very  thin  fingers  of  tenuous  gas  mov- 
ing into  the  denser  region.  The  separation  between  dense  columns  decreases  as  the 
SGI  evolves.  The  growth  rates  for  RT  models  are  determined  with  the  same  In- 
linear  plotting  method  used  for  the  SGI  cases.  However,  for  RT  models  an  upward 
moving  dense  column  is  used.  The  growth  rate  for  model  hart257dl  (figure  3.16)  is 
a;  = 0.580tj^.  This  result  is  surprisingly  low  considering  how  well  the  experimental 
growth  rate  matches  the  theoretical  growth  rate  for  the  SGI  models  so  far  considered. 
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Figure  3.14:  Velocity  growth  for  model  hnsg257dl.  Data  points  are  normalized  to 
the  velocity  value  at  to  = l.Otg-  The  best-fit  line  using  data  points  from  O.SOtg  to 
l.SOte  yields  a growth  rate  of  1.031^7^  The  e-folding  time  is  tg  = 8.458  x 10^^  s. 
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Figure  3.15:  Evolution  of  the  RT  instability  for  model  hart257dl.  Filled  contours 
show  the  density  distribution  at  a)lte;  b)2te;  c)3te;  and  d)4te-  The  e-folding  time  is 
te  = 8.458  X 10^^  s.  The  cell  size  Ax  = Ay  = 1.0  x 10^®cm. 
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Figure  3.16:  Velocity  growth  for  RT  model  hart257dl.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l-Otg-  The  best-fit  line  using  data  points  from  O.SOtg  to 
l.SOte  yields  a growth  rate  of  0.580t~^.  The  e-folding  time  is  — 8.458  x 10^^  s. 
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Models  for  both  types  of  instabilities  deviate  from  the  theoretical  assumptions  of  in- 
compressibility and  non-stratification.  The  lower  than  expected  RT  growth  rate  may 
be  attributable  to  the  use  of  2-D  simulations  to  model  a 3-D  phenomenon.  Calder  et 
al.  (2000)  have  demonstrated  that  RT  simulations  performed  with  3-D  perturbations 
produce  significantly  higher  growth  rates  than  those  with  2-D  perturbations.  Three 
dimensional  simulations  of  the  SGI  are  needed  to  verify  that  the  same  is  not  true  in 
the  self-gravitating  case. 

According  to  equation  (3.11),  a higher  density  contrast  at  the  interface  results 
in  a faster  growth  rate  for  both  the  SGI  and  the  RT  instability.  Models  hnsg257d2  and 
hart257d2  are  high  temperature  models  with  pi/p2  = 5.  Region  1 remains  unchanged 
from  models  hnsg257dl  and  hart257dl,  but  region  2 is  assigned  the  interfacial  values 
P2  = 4 X 10~^^p  cm~^  and  T2  = 5000  K.  The  perturbation  velocity  and  J value  are 
defined  using  values  in  region  1,  and  hence  remain  unchanged.  The  higher  density 
contrast  results  in  a shorter  theoretical  e-folding  time,  tg  = 4.728  x lO^^s.  Figures 
(3.17)  and  (3.18)  show  the  growth  of  the  instabilities  for  model  hnsg257d2  and  model 
hart257d2,  respectively.  Note  that  the  panels  show  the  density  distribution  at  dif- 
ferent times  than  in  the  previous  figures.  The  larger  density  contrast  at  the  interface 
results  in  a larger  momentum  difference,  so  penetration  by  the  tenuous  fingers  into 
the  denser  medium  is  impeded  early  in  the  simulation.  Compared  to  the  pi/ p2  — 2 
models,  both  the  SGI  and  the  RT  instability  grow  faster  in  terms  of  absolute  time.  In 
terms  of  the  theoretical  e-folding  time,  the  SGI  growth  is  slower.  The  growth  rates  for 
the  SGI  and  RT  models  are  far  more  similar  than  in  the  previous  case.  The  velocity 
growth  for  model  hnsg257d2  (figure  3.19)  yields  a growth  rate  of  u = 0.735t~^.  The 
velocity  growth  for  model  hart257d2  (figure  3.20)  yields  a growth  rate  of  u;  = 0.720tj^. 
Both  models  have  growth  rates  lower  than  the  theoretical  prediction. 

The  density  ratio  of  5 was  chosen  because  it  produces  an  e-folding  time  com- 
parable to  the  free-fall  time  of  the  denser  medium.  A sphere  with  the  density  and 
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Figure  3.17:  Evolution  of  the  SGI  for  model  hnsg257d2.  Filled  contours  show  the 
density  distribution  at  a)2tg;  b)3te;  c)4te;  and  d)5te-  The  e- folding  time  is  tg  = 
4.728  X 10^^  s.  The  cell  size  Ax^Ay  — 1.0  x 10^®cm. 
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Figure  3.18:  Evolution  of  the  RT  instability  for  model  hart257d2.  Filled  contours 
show  the  density  distribution  at  a)2te]  b)3tg;  c)4te;  and  d)5te-  The  e-folding  time  is 
te  = 4.728  X 10^^  s.  The  cell  size  Ax  = Ay  = 1.0  x lO^^cm. 
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Figure  3.19:  Velocity  growth  for  SGI  model  hnsg257d2.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l-Ote.  The  best-fit  line  using  data  points  from  0.70te  to 
2.00te  yields  a growth  rate  of  0.735^7^  The  e-folding  time  is  te  — 4.728  x 10^^  s. 
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Figure  3.20:  Velocity  growth  for  RT  model  hart257d2.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l-0te-  The  best-fit  line  using  data  points  from  0.70te  to 
2.00te  yields  a growth  rate  of  0.720t~h  The  e-folding  time  is  tg  = 4.728  x 10^^  s. 
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temperature  of  region  1 in  the  high  temperature  models  has  a free-fall  time  (equation 
2.41)  of  = 4.697  X lO^^s.  In  the  ISM,  density  contrasts  at  the  boundaries  of  clouds 
and  filaments  are  typically  far  greater  than  5.  In  these  regions,  the  SGI  can  grow  on 
time-scales  far  shorter  than  the  free-fall  time. 

3.4.3  Low  Temperature  Models 

The  low  temperature  models  have  density  and  temperature  values  appropriate 
for  cold  molecular  clouds,  p\  = and  T\  = 10  K.  The  higher  density  and 

lower  temperature  result  in  lower  values  for  the  scale  height  and  Jeans  length  than 
in  the  hot  models.  Therefore,  if  given  the  same  perturbation  wavelength  as  the 
hot  models,  the  cold  models  have  higher  J values.  Whereas  the  hot  models  have 
a J value  of  0.0174,  some  of  the  cold  models  have  J values  as  high  as  0.388.  As 
a result,  the  cold  models  begin  with  a more  stratified  density  distribution  and  are 
more  likely  to  undergo  compression.  The  sound  speed  is  lower  for  the  cold  models 
(ci  = 24, 124  cm  s“^)  so  a larger  fraction  of  the  sound  speed  Uq/ci  = 0.05  is  used  as 
the  velocity  perturbation  amplitude.  A greater  departure  from  the  incompressible, 
non-stratified  theoretical  growth  rate  is  expected. 

The  first  set  of  cold  models  has  a density  ratio  of  5 and  a high  J value.  The 
models  with  the  highest  J value  (J  = 0.388)  are  the  SGI  model  cnsg257el  and  the  RT 
model  cart257el.  As  the  most  stratified  models,  they  begin  with  a larger  departure 
from  static  equilibrium  than  any  other  model.  Even  so,  the  motions  induced  by  the 
departure  from  equilibrium  are  not  large  enough  to  affect  the  instability  growth  at 
the  interface.  The  theoretical  e-folding  time  for  these  models  is  te  = 2.115  x lO^^s.  The 
growth  rate  is  higher  than  in  previous  cases  due  to  the  higher  background  density. 
The  evolution  of  model  cnsg257el  is  shown  in  figure  (3.21).  Only  two  full  cosine 
waves  are  excited.  As  previously  seen,  the  tenuous  gas  streams  into  the  denser  gas 
through  narrow  channels,  leaving  wide,  dense  columns  behind.  By  the  tenuous 
columns  have  hit  the  bottom  boundary  of  the  grid.  The  RT  model  cart257el  is 
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Figure  3.21;  Evolution  of  the  SGI  for  model  cnsg257el.  Filled  contours  show  the 
density  distribution  at  a)3tg;  b)4tg;  c)5te;  and  d)6te-  Note  the  resemblance  to  the  RT 
model  cart257el  at  late  times.  The  e- folding  time  is  tg  = 2.115  x 10^^  s.  The  cell  size 
Ax  = Ay  = 0.5  X lO^^cm. 
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Figure  3.22:  Evolution  of  the  RT  instability  for  model  cart257el.  Filled  contours 
show  the  density  distribution  at  a)3te;  b)4te;  c)5te]  and  d)6te-  Note  the  resemblance 
to  the  SGI  model  cnsg257el  at  late  times.  The  e-folding  time  is  tg  = 2.115  x 10^^  s. 
The  cell  size  Aa:  = Ay=:0.5  x 10^®cm. 
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shown  in  figure  (3.22).  As  expected,  more  compression  is  seen  in  the  cold  SGI  and 
RT  models  than  in  the  hot  models.  As  in  other  RT  models,  the  dominant  features  in 
model  cart257el  are  dense  columns  pushing  up  and  displacing  the  tenuous  fluid.  Note 
the  striking  similarities  between  the  morphologies  of  models  cnsg257el  and  cart257el. 
Although  the  dynamic  processes  are  different  and  the  instabilities  grow  in  different 
ways,  both  give  rise  to  dense  columns  of  comparable  size  and  width.  In  both  cases,  the 
height  of  the  dense  columns  after  6te  is  roughly  170  cells  or  8.5  x 10^^cm  = 0.275pc. 
The  largest  difference  is  the  size  of  the  mushroom  cap  at  the  head  of  each  dense 
column.  Whichever  fluid  moves  fastest,  the  tenuous  for  SGI  and  the  dense  for  RT, 
forms  narrower  columns  and  smaller  caps  than  the  slower  moving  fluid.  As  in  the 
high  temperature  models  with  pi/p2  = 5,  the  growth  rates  for  models  cnsg257el  and 
cart257el  are  in  close  agreement.  The  velocity  growth  for  model  cnsg257el  (figure 
3.23)  yields  a growth  rate  of  a;  = 0.852t“^.  The  velocity  growth  for  model  cart257el 
(figure  3.24)  yields  a growth  rate  of  o;  = 0.836t7^.  Unexpectedly,  the  experimental 
growth  rate  for  this  highly  stratified  RT  model  shows  better  agreement  with  the  non- 
stratified  theoretical  prediction  than  the  growth  rates  for  the  high  temperature  RT 
models. 

The  prediction  that  the  SGI  growth  rate  is  independent  from  the  perturba- 
tion wavelength  was  put  through  another  test  with  low  temperature  models.  Lower 
resolution  models  were  excited  with  2,  3,  4,  and  5 waves.  The  wave  numbers  and  J 
values  varied  accordingly.  The  result  is  shown  in  figure  (3.25).  As  in  the  moderate 
temperature  version,  ky  remained  constant  amongst  the  models  while  k^  changed. 
The  models  with  the  longer  perturbation  wavelengths  have  slightly  higher  growth 
rates,  not  a surprising  result  considering  the  range  of  J values  represented. 

Additional  simulations  were  performed  using  low  temperature  models  with  a 
density  ratio  of  5 to  explore  the  effects  of  varying  the  J value.  The  physical  parame- 
ters are  unaltered,  but  a smaller  cell  size  is  used.  For  example,  model  cnsg257d2  has 
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Figure  3.23:  Velocity  growth  for  SGI  model  cnsg257el.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l.Ote-  The  best-fit  line  using  data  points  from  l.Ote  to 
S.Ote  yields  a growth  rate  of  0.852t“F  The  e-folding  time  is  tg  = 2.115  x 10^^  s. 
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Figure  3.24;  Velocity  growth  for  RT  model  cart257el.  Data  points  are  normalized  to 
the  velocity  value  at  tQ  = lMe-  The  best-fit  line  using  data  points  from  l.Otg  to  S.Ote 
yields  a growth  rate  of  O.SSGt^h  The  e-folding  time  is  tg  — 2.115  x 10^^  s. 


77 


60 

40 

20 

0 

-20 

-40 

-60 


cnsgl 29e4 


density  (x10‘^’) 
o = 10.0 
b=7.21 
c=4.42 
d = 1.62 


I ...  I ■_  ■.■.■■L,.  ■■  I ...  I ...  I ...  I 


-60  -40  -20  0 20  40  60 

(d) 


Figure  3.25:  Comparison  of  cold  SGI  models  with  varying  k^.  The  panels  show  the 
density  contours  at  t = 4te  for  models  a)cnsgl29el;  b)cnsgl29e2;  c)cnsgl29e3;  and 
d)cnsgl29e4.  Longer  wavelength  models  show  a slightly  larger  growth  rate.  The 
e-folding  time  is  te  = 2.115  x 10^^  s for  all  four  models.  The  cell  size  Aa;  = Ay  = 
1.0  X 10^®cm. 
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a cell  size  one-fifth  as  large  as  the  cell  size  in  model  cnsg257el.  The  perturbation 
wavelength  and  J value  (J  = 0.0776)  are  also  one-fifth  as  large.  The  result  is  shown 
in  figure  (3.26).  The  same  velocity  perturbation  is  used  for  models  cnsg257d2  and 
cnsg257el.  The  fluid  travels  across  the  grid  in  less  time  due  to  the  smaller  cell  size. 
After  a time  of  3.5te,  the  dense  columns  are  about  140  cells  tall,  which  is  a physical 
distance  of  only  1.4x  10^^cm  = 0.0454pc.  By  4.5te  the  tenuous  fingers  have  rebounded 
from  the  bottom  boundary  of  the  grid.  For  the  velocities  used,  the  physical  extent  of 
the  grid  is  not  large  enough  to  follow  the  SGI  growth  for  very  long.  Considering  this 
restriction,  model  cnsg257d2  does  a fair  job  of  reproducing  the  interfacial  region  of 
model  cnsg257el.  The  growth  rate  for  model  cnsg257d2  (figure  3.27)  is  a;  = 0.787tj^, 
a value  close  to  that  for  model  cnsg257el.  (The  velocity  growth  in  model  cnsg257d2 
does  not  result  in  a clearly  defined  best-fit  line  on  the  In-linear  growth  plot.  The  ex- 
perimental value  of  the  growth  rate  is  more  uncertain  than  that  of  model  cnsg257el.) 
Despite  the  differences  in  cell  size  and  J value,  models  cnsg257el  and  cnsg257d2  pro- 
duce comparable  results.  Contrary  to  expectations,  a larger  degree  of  stratification 
does  not  result  in  a greater  deviation  from  the  theoretical  growth  rate. 

In  the  low  temperature,  high  J regime  (models  cnsg257el  and  cart257el), 
the  SGI  and  the  RT  instability  resulted  in  very  similar  morphologies  and  growth 
rates.  Such  a similarity  is  not  readily  seen  in  the  low  temperature,  low  J models. 
The  low  J model  cart257d2  is  shown  in  figure  (3.28).  The  dense  columns  in  model 
cart257d2  are  of  comparable  height  to  those  in  model  cnsg257d2  throughout  most  of 
the  simulation.  However,  there  is  a greater  discrepancy  in  the  size  of  the  mushroom 
caps.  The  constant  acceleration  g in  model  cart257d2  is  one-fifth  that  in  model 
cart257el.  The  two  cold  RT  models  show  very  similar  structures,  but  the  growth 
rates  are  quite  different.  The  velocity  growth  for  model  cart257d2  (figure  3.29)  yields 
a growth  rate  co  = 0.465tg,  much  lower  than  the  growth  rate  for  model  cart257el. 
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Figure  3.26:  Evolution  of  the  SGI  for  model  cnsg257d2.  Filled  contours  show  the 
density  distribution  at  a)1.5te;  b)2.5te;  c)3.5tg;  and  d)4.5tg.  The  e-folding  time  is 
tg  = 2.115  X 10^^  s.  The  cell  size  Aa:  = Ay  = 0.1  x lO^^cm. 
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Figure  3.27:  Velocity  growth  for  SGI  model  cnsg257d2.  Data  points  are  normalized 
to  the  velocity  value  at  to  = 1-Ote-  The  best-fit  line  using  data  points  from  1.2te  to 
2.6tg  yields  a growth  rate  of  0.787^7^  The  e-folding  time  is  te  = 2.115  x 10^^  s. 
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(c)  (d) 

Figure  3.28:  Evolution  of  the  RT  instability  for  model  cart257d2.  Filled  contours 
show  the  density  distribution  at  a)1.5te;  b)2.5te;  c)3.5te]  and  d)4.5tg-  The  e-folding 
time  is  tg  = 2.115  x 10^^  s.  The  cell  size  Ax  = Ay  = 0.1  x 10^®cm. 
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Figure  3.29:  Velocity  growth  for  RT  model  cart257d2.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l.Ote-  The  best-fit  line  using  data  points  from  l.Ote  to 
3.0te  yields  a growth  rate  of  0.465tjh  The  e-folding  time  is  tg  = 2.115  x 10^^  s. 
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Model  cart257d2  resembles  the  high  temperature  pi/p2  = 2 models  in  that  its  growth 
rate  is  lower  than  the  corresponding  SGI  model  cnsg257d2. 

The  low  temperature,  low  J models  were  also  performed  with  a perturbation 
that  excited  four  full  cosine  waves.  Comparison  with  the  previous  set  of  models 
provides  some  idea  of  how  the  spacing  of  the  peaks  on  the  grid  affects  the  evolution 
of  the  instabilities.  The  SGI  model  cnsg257dl,  shown  in  figure  (3.30),  is  similar  to 
model  cnsg257d2,  except  that  the  values  of  kx  and  ky  are  twice  as  big.  In  addition 
to  reducing  J by  one-half,  the  increase  in  k results  in  a steeper  exponential  falloff 
in  velocity  away  from  the  interface,  so  the  initial  perturbation  amplitude  is  less. 
This  difference  notwithstanding,  the  morphologies  are  not  very  different  until  late  in 
the  evolution.  By  4tg,  model  cnsg257dl  shows  signs  of  asymmetries  in  the  tenuous 
columns.  The  fingers  of  gas  have  begun  interacting.  The  cell  size  is  too  small  to  fit 
four  non-interacting  waves  across  the  grid.  When  only  two  waves  are  perturbed,  the 
columns  are  separated  by  a sufficient  distance  so  that  interaction  between  them  is 
delayed.  The  growth  rate  for  model  cnsg257dl  (figure  3.31)  is  a bit  problematic  to 
calculate  due  to  the  interaction  of  the  gas  streams.  Limiting  the  slope  calculation  to 
data  points  between  Itg  and  2te  results  in  a growth  rate  of  cj  = 0.620te,  a value  lower 
than  that  for  model  cnsg257d2. 

The  RT  model  cart257dl  with  four  excited  waves  compares  to  model  cart257d2 
in  very  much  the  same  way.  Because  the  RT  growth  rate  depends  upon  the  perturba- 
tion wavelength,  a smaller  value  of  g was  used  for  model  cart257dl.  Figures  (3.32)  and 
(3.33)  show  the  morphology  and  velocity  growth,  respectively,  for  model  cart257dl. 
The  morphologies  of  model  cart257dl  and  cart257d2  are  similar,  but  the  growth  rate 
for  model  cart257dl  is  slightly  smaller  than  that  for  model  cart257d2.  Unlike  the  SGI 
model,  no  interaction  between  the  growing  dense  columns  is  seen  in  model  cart257dl. 
The  reason  is  partly  due  to  the  smaller  growth  rate  for  model  cart257dl,  compared  to 
that  for  model  cnsg257dl,  and  partly  due  to  the  lack  of  self-gravity  induced  motion 
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Figure  3.30:  Evolution  of  the  SGI  for  model  cnsg257dl.  Filled  contours  show  the 
density  distribution  at  a)2tg;  b)3tg;  c)4te;  and  d)5te-  The  e-folding  time  is  te  = 
2.115  X 10^^  s.  The  cell  size  Aa:  = Ar/  = 0.1  x 10^®cm. 
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Figure  3.31;  Velocity  growth  for  SGI  model  cnsg257dl.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l.Otg-  The  best-fit  line  using  data  points  from  l.Ote  to 
2.0te  yields  a growth  rate  of  0.620^7^  The  e-folding  time  is  te  = 2.115  x 10^^  s. 
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in  the  x direction.  For  both  the  SGI  and  RT  models  considered,  altering  the  num- 
ber of  perturbed  waves  across  the  grid  has  little  effect,  provided  the  spacing  is  large 
enough  (and  small  enough)  to  prevent  interaction  between  the  columns.  Having 
fewer  waves  on  the  grid  means  more  cells  make  up  the  width  of  each  column  and  cap, 
which  allows  more  detailed  structures  to  be  seen. 

All  of  the  low  temperature  models  so  far  mentioned  have  a density  ratio  of  5 
at  the  interface.  Model  cnsg257d3  has  a density  ratio  of  2 at  the  interface,  resulting 
in  a longer  e-folding  time.  The  growth  of  the  SGI  is  shown  in  figure  (3.34).  The 
lower  momentum  contrast  at  the  interface  allows  the  SGI  to  grow  sooner  relative  to 
the  e-folding  time  than  when  pi/p2  = 5.  The  tenuous  gas  streams  strike  the  bottom 
of  the  grid  rather  quickly.  The  structure  seen  in  model  cnsg257d3  resembles  that  in 
the  high  temperature  models  with  pxjp^  — l.  The  tenuous  fingers  become  narrower 
as  the  simulation  progresses.  The  structure  at  later  times  is  very  interesting.  Figure 
(3.35)  shows  model  cnsg257d3  at  a time  of  5tg.  Secondary  instabilities  have  formed 
both  on  the  dense  columns  and  between  them.  The  growth  rate  for  this  model 
(figure  3.36)  is  a;  = 0.882,  which  is  closer  to  the  theoretical  prediction  than  any  other 
low  temperature  model.  However,  the  growth  rate  for  the  corresponding  RT  model 
cart257d3  (not  shown)  is  again  unexpectedly  low  at  a;  = 0.306.  The  hot  and  cold 
Pi/P2  = 2 models  agree  with  one  another  fairly  well.  The  hot  models  were  not  run 
long  enough  to  see  any  structures  akin  to  figure  (3.35),  but  this  does  not  preclude 
the  possibility  such  morphologies  do  arise  at  late  times.  The  morphology  of  model 
cnsg257d3  at  t = 3.5te  is  not  very  different  from  the  morphology  of  model  hnsg257dl 
at  t = 5.0te  (figure  3.13).  The  largest  difference  between  the  hot  and  cold  models  is 
the  time-scale  on  which  they  grow,  but  the  difference  in  e-folding  time  is  due  to  the 
change  in  density,  not  the  change  in  temperature. 
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Figure  3.32:  Evolution  of  the  RT  instability  for  model  cart257dl.  Filled  contours 
show  the  density  distribution  at  a)2tg;  h)3te]  c)4tg;  and  d)5te-  The  e-folding  time  is 
te  = 2.115  X 10^^  s.  The  cell  size  Ax  = Ay  = 0.1  x 10^®cm. 


88 


0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4  2.6  2.8  3.0 

time  in  e— folding  times 


Figure  3.33:  Velocity  growth  for  RT  model  cart257dl.  Data  points  are  normalized 
to  the  velocity  value  at  to  = 1-Otg-  The  best-fit  line  using  data  points  from  l.Otg  to 
2.0tg  yields  a growth  rate  of  0.397tJ^.  The  e-folding  time  is  tg  — 2.115  x 10^^  s. 
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Figure  3.34:  Evolution  of  the  SGI  for  model  cnsg257d3.  Filled  contours  show  the 
density  distribution  at  a)0.5te;  b)1.5te;  c)2.5te;  and  d)3.5te.  The  e-folding  time  is 
te  = 3.783  X 10^^  s.  The  cell  size  Ax  = Ay  = 0.1  x 10^®cm. 


Y cell 


90 


-100  -50  0 50  100 

X cell 


Figure  3.35:  Model  cnsg257d3  at  t = hte.  The  tenuous  fingers  have  long  since  re- 
bounded from  the  bottom  boundary,  but  the  dense  gas  is  only  beginning  to  impact 
the  top  boundary. 
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Figure  3.36:  Velocity  growth  for  SGI  model  cnsg257d3.  Data  points  are  normalized 
to  the  velocity  value  at  to  = l.Otg.  The  best-fit  line  using  data  points  from  0.7tg  to 
1.6tg  yields  a growth  rate  of  0.882^7^.  The  e-folding  time  is  tg  — 3.783  x 10^^  s. 
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3.5  Discussion 

The  work  presented  in  this  chapter  represents  the  first  computational  study 
of  the  SGI.  Although  the  parameter  space  considered  is  not  very  large  (in  particular, 
models  with  higher  density  ratios  are  desirable)  some  general  remarks  can  be  made. 
The  SGI  is  in  some  ways  similar  to  the  RT  instability,  but  in  other  ways  the  two  are 
very  different.  Both  destabilize  an  interface  of  discontinuous  density;  the  SGI  in  the 
presence  of  self-gravity,  and  the  RT  instability  when  a dense  gas  is  accelerated  into 
a tenuous  gas.  Both  instabilities  preserve  the  perturbation  wavelength  well  into  the 
nonlinear  regime.  However,  the  two  instabilities  grow  in  very  different  manners.  The 
growth  of  the  SGI  is  characterized  by  thin  fingers  of  tenuous  material  streaming  into 
the  denser  material.  In  the  case  of  an  interface  with  a low  density  ratio,  the  dense 
medium  is  displaced  very  little.  The  RT  instability  develops  in  the  opposite  manner; 
thin,  dense  columns  stream  into  the  tenuous  medium.  When  the  RT  instability 
acts  upon  interfaces  with  different  density  ratios,  the  growth  rate  changes,  but  the 
resultant  morphologies  of  the  evolved  gas  remain  very  much  the  same.  The  SGI  gives 
rise  to  very  different  morphologies  when  the  density  ratio  is  altered.  This  will  be 
explored  further  in  future  work.  A careful  look  at  the  velocity  growth  plots  reveals 
another  difference  between  the  SGI  and  the  RT  instability.  After  an  early  interval 
of  exponential  growth,  the  RT  growth  rate  decreases  in  the  nonlinear  stage  (e.g. 
Youngs  1984).  Conversely,  the  growth  of  the  SGI  accelerates  in  the  nonlinear  phase. 
(This  conclusion  must  be  tempered  by  the  reminder  that  the  boundary  conditions 
used  by  the  gravity  solver  become  less  appropriate  as  the  instability  develops.)  For 
all  of  their  differences,  the  SGI  and  the  RT  instability  can  produce  structures  with 
similar  morphologies  when  the  density  ratio  at  the  interface  is  large.  This  leads  to 
the  speculation  that  dense  columns  of  gas  can  form  in  the  ISM  due  to  either  the  SGI 
or  RT  instability.  The  surrounding  environment  and  the  nature  of  the  perturbation 
determine  which  instability  is  dominant. 
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Table  3.3:  Experimental  growth  rates  for  SGI  and  RT  models 


model 

Pl/P2 

UJ 

tl 

hasg257dl 

2 

1.028tji 

8.458 

hnsg257dl 

2 

1.03R-1 

8.458 

hnsg257d2 

5 

0.735f-^ 

4.728 

hcirt257dl 

2 

0.580tji 

8.458 

hart257d2 

5 

0.720t-^ 

4.728 

cnsg257el 

5 

0.852tji 

2.115 

cnsg257dl 

5 

0.620t“i 

2.115 

cnsg257d2 

5 

0.787^7^ 

2.115 

cnsg257d3 

2 

0.882^^1 

3.783 

cart257el 

5 

0.836t7^ 

2.115 

cart257dl 

5 

0.397^7! 

2.115 

cart257d2 

5 

0.465t"^ 

2.115 

cart257d3 

2 

0.306tji 

3.783 

is  given  in  units  of  10^ 


The  theory  of  HWL  does  a very  good  job  of  predicting  growth  rates  for  the 
SGI.  Even  though  the  growth  rates  given  by  equation  (3.11)  were  derived  in  the 
incompressible,  non-stratified  limit,  the  experimental  growth  rates  agree  fairly  well. 
The  experimental  growth  rates  are  summarized  in  Table  (3.3).  For  the  SGI  models, 
the  experimental  growth  rates  vary  from  0.620tg  to  1.031te-  The  SGI  models  with 
an  interfacial  density  ratio  of  2 consistently  come  closer  to  matching  the  theoretical 
growth  rate  than  models  with  a density  ratio  of  5.  The  experimental  growth  rates  for 
the  low  temperature  SGI  models  do  not  show  a strong  dependence  on  the  J value, 
although  more  examples  are  needed  to  confirm  this.  The  same  cannot  be  said  of  the 
RT  models.  The  experimental  RT  growth  rates  are  all  lower  than  the  theoretical 
prediction.  Contrary  to  expectations,  the  closest  agreement  between  theoretical  and 
experimental  values  is  seen  for  the  RT  model  with  the  largest  degree  of  stratification. 
In  all  cases,  the  RT  growth  rate  is  lower  than  the  corresponding  SGI  growth  rate, 
perhaps  indicating  that  the  SGI  growth  is  impacted  less  by  the  use  of  2-D  simulations. 

Very  little  difference  is  seen  between  the  7 = 2 and  the  7 = 1.4  models.  This 
independence  from  7 may  not  hold  up  as  7 approaches  unity  or  the  velocity  amplitude 
approaches  the  sound  speed.  The  grid  resolution  has  little  effect  upon  the  growth 
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rates  of  the  instabilities,  but  a higher  resolution  allows  more  detailed  structures  to  be 
seen.  The  best  cell  size  for  problems  of  this  sort  is  one  that  strikes  a balance  between 
the  need  for  adequate  spatial  resolution  and  the  desire  to  minimize  boundary  effects. 

Far  more  work  is  needed  to  understand  the  SGI  and  its  importance  in  the 
ISM.  Many  variations  of  hydrostatic  models  are  possible.  Different  values  of  density, 
density  ratio,  temperature,  7,  J,  and  velocity  perturbation  amplitude  remain  to 
be  explored.  This  study  is  limited  to  planar  symmetry  in  Cartesian  coordinates. 
Given  the  plethora  of  filamentary  structures  seen  in  the  ISM,  the  development  of 
the  SGI  in  cylindrical  geometry  bears  investigation.  The  SGI  behaves  differently  in 
cylindrical  geometry  (Hunter  et  al.  1998).  A cylindrical  interface  at  a radial  distance 
To  with  interior  density  p\  and  exterior  density  p2  is  stable  to  kink  modes  (azimuthal 
perturbations)  but  unstable  to  sausage  modes  (radial  perturbations).  Unlike  the 
planar  case,  growth  rates  for  the  sausage  modes  depend  upon  the  wave  number  of 
the  perturbation.  There  exists  a maximum  wave  number,  or  minimum  perturbation 
wavelength,  beyond  which  a cylindrical  interface  is  stable  to  the  SGI.  The  condition 
for  instability  is  kro  < 1.17(1—  p2/pi)-  Future  work  will  include  the  modeling  of  the 
SGI  along  cylindrical  interfaces. 

As  stated  earlier  in  this  chapter,  assuming  an  initial  state  of  hydrostatic  equi- 
librium is  not  appropriate  for  the  ISM.  A thorough  study  of  the  ISM  must  include 
magnetic  fields,  turbulence,  velocity  streams,  shocks,  heating,  and  cooling  along  with 
gravitational  forces.  Each  process  alters  the  morphology  and  evolution  of  the  ISM, 
often  in  different  ways  and  on  different  time-scales.  As  the  understanding  of  the  SGI 
improves,  additional  physics  will  be  added  to  the  models  in  order  to  form  a more 
accurate  picture  of  the  ISM.  The  relative  importance  of  the  SGI,  when  coupled  with 
other  instabilities,  remains  to  be  determined. 


CHAPTER  4 

COLLIDING  GAS  FLOWS  IN  GARTESIAN  GEOMETRY 


The  ISM  is  a dynamic  medium  permeated  by  supersonic  flows,  which  give  rise 
to  shocks.  Since  shocks  are  always  compressive,  density  enhancements  result.  This 
has  long  fueled  speculation  that  converging  velocity  fields  can  trigger  protostellar 
collapse  in  regions  that  do  not  initially  contain  a Jeans  mass  of  material  (e.g.  Hunter 
1979,  Hunter  and  Fleck  1982,  Tohline  et  al.  1987).  In  addition  to  their  possible  role  in 
facilitating  the  star  formation  process,  interstellar  shocks  are  worth  investigating  for 
their  effects  on  cloud  morphology.  Shocks  are  prone  to  various  instabilities  (Draine 
and  McKee  1993)  which  dramatically  alter  the  shape  of  the  ISM.  In  this  chapter  and 
the  next,  I consider  the  instability  of  a dense  slab  bounded  by  accretion  shocks,  a 
problem  previously  studied  by  Hunter,  Sandford,  Whitaker,  and  Klein  (1986,  here- 
after HSWK).  At  issue  is  whether  the  head-on  collision  of  supersonic  flows  can  lead 
to  structures  which  fragment  into  self-gravitating  clumps.  Simulations  of  colliding 
flows  were  performed  in  two  dimensions  using  both  Gartesian  and  cylindrical  grid- 
s.  Each  has  its  advantages  and  disadvantages.  A cylindrical  grid  is  better  suited  for 
modeling  filamentary  structures  or  systems  in  which  gas  flows  converge  at  oblique  an- 
gles. Self-gravity  pulls  matter  towards  the  symmetry  axis,  forming  structures  which 
approximate  three  dimensional  spheres.  Thus,  the  most  interesting  structures  (self- 
gravitating  clumps)  form  along  the  grid  boundary  representing  the  symmetry  axis. 
The  boundary  conditions  must  be  carefully  controlled.  Models  in  Gartesian  geometry 
have  no  preferred  direction.  However,  a Cartesian  grid  cannot  truly  represent  any 
3-D  structure  since  the  third  dimension  is  infinite.  Self-gravitating  clumps  are  akin 
to  infinitely  extended,  radially  collapsing  cylinders. 
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Cylindrical  models  are  presented  in  the  next  chapter.  In  this  chapter,  Carte- 
sian models  are  considered.  First,  I summarize  the  studies  of  colliding  flows  and 
clouds  previously  done  by  others.  Then,  the  Cartesian  models  are  described  and 
results  presented.  The  chapter  ends  with  a discussion  of  how  the  current  results 
compare  to  previous  studies. 

4.1  Background 

Astronomers  have  long  known  that  any  accurate  model  of  the  ISM  must  ac- 
count for  its  dynamical  nature  (e.g.  von  Weizsacker  1951,  Chandrasekhar  and  Miinch 
1952).  As  numerical  methods  have  improved  over  time,  many  attempts  have  been 
made  to  model  dynamical  processes  in  the  ISM.  A variety  of  studies  have  considered 
the  collision  of  two  supersonic  gas  streams,  the  collision  of  two  clouds,  or  the  interac- 
tion between  a stream  and  a cloud.  In  all  cases,  collisions  give  rise  to  dense  slabs  or 
disks  bounded  by  accretion  shocks.  At  question  is  whether  the  interactions  of  streams 
or  clouds  tends  to  enhance  or  inhibit  star  formation.  In  an  early  computational  study 
of  colliding  clouds.  Stone  (1970a,  1970b)  found  that  although  the  collision  resulted  in 
a large  compression,  the  clouds  ultimately  re-expanded  before  gravitational  collapse 
could  occur.  Hunter  and  Fleck  (1982)  reached  a different  conclusion.  They  used  a 
virial  approach  to  show  that  interacting  flows  enhanced  star  formation.  Regions  that 
start  off  with  an  insufficient  mass  for  collapse  can  be  compressed  by  convergent  flows 
with  efficient  cooling.  The  higher  density  leads  to  a smaller  Jeans  mass,  allowing 
collapse  to  commence.  HSWK  followed  up  on  this  idea  with  computational  models 
of  colliding  supersonic  gas  streams.  The  models,  performed  in  cylindrical  geometry, 
included  self-gravity  and  molecular  cooling.  The  colliding  flows  gave  rise  to  a dense 
disk  bounded  by  accretion  shocks.  The  disk  was  found  to  be  unstable  due  to  the 
combined  effects  of  efficient  cooling  and  the  continuous  ram  pressure  provided  by  the 
accreting  shocks.  The  instability  caused  the  disk  to  fragment  into  dense  clumps.  In 
the  absence  of  self-gravity,  the  clumps  were  short  lived  objects.  When  gravity  was 
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included  in  the  simulation,  fragments  along  the  symmetry  axis  continued  to  accrete 
gas  until  surpassing  the  local  Jeans  mass.  These  results  support  the  hypothesis  that 
colliding  streams  enhance  the  star  formation  process. 

The  instability  discovered  by  HSWK  was  subsequently  observed  in  simulations 
of  colliding  wind  binaries  by  Stevens,  Blondin,  and  Pollack  (1992,  hereafter  SBP), 
who  described  it  as  a thin  shell  instability.  Although  they  considered  much  hotter 
gas  and  higher  flow  speeds  than  HSWK,  SBP  found  the  morphology  of  the  thin 
shell  instability  to  be  very  similar  to  the  earlier  models.  They  also  reiterated  the 
need  for  a cooling  mechanism  in  the  post-shock  slab.  The  computational  results 
prompted  Vishniac  (1994)  to  perform  an  analytic  study  of  the  instability.  He  found 
that  when  isothermal  flows  collide,  the  slab  is  stable  to  linear  order,  but  unstable  to 
nonlinear  effects  once  the  displacement  is  comparable  to  the  thickness  of  the  slab. 
He  termed  the  instability  the  “nonlinear  thin  shell  instability”  (NTSI).  The  physical 
basis  for  the  NTSI  is  the  transfer  of  momentum  such  that  perturbations  are  enhanced. 
Suppose  that  a dense  slab  is  displaced  by  a bending  mode  perturbation.  Gas  which 
strikes  the  bottom  of  the  slab  will  be  transported  transversely  into  the  regions  which 
are  displaced  upward.  Likewise,  the  gas  striking  the  top  of  the  perturbed  slab  will 
transfer  some  of  its  momentum  to  the  downward  moving  peaks.  The  cumulative  effect 
causes  the  displacement  to  grow.  The  momentum  transfer  fueling  the  growth  of  the 
instability  must  compete  against  the  stabilizing  effects  of  enhanced  ram  pressure  on 
the  peaks  moving  into  the  accreting  flows  and  the  tendency  of  the  post-shock  flow  to 
All  in  the  ripples.  Linear  instability  is  suppressed,  so  only  the  nonlinear  instability 
(finite  initial  amplitude)  can  grow.  Vishniac  (1994)  derived  an  approximate  growth 
rate  for  the  NTSI 

u^Csk(kL)^,  (4.1) 

where  is  the  sound  speed  in  the  slab,  k the  perturbation  wavenumber,  and  L the 
slab  displacement.  Vishniac  (1994)  confined  his  analysis  to  isothermal  flows,  and  was 


98 


therefore  unable  to  comment  on  the  consequences  of  altering  cooling  rates.  He  did 
however  indicate  that  a high  compression  is  needed  for  the  development  of  the  NTSI. 
An  efficient  cooling  mechanism  gives  rise  to  higher  compressions  than  the  isothermal 
case,  and  could  therefore  enhance  the  NTSI.  More  compression  will  also  result  from 
more  highly  supersonic  flows. 

The  most  comprehensive  computational  study  of  the  NTSI  was  performed 
by  Klein  and  Woods  (1998,  hereafter  KW).  They  used  adaptive  mesh  refinement  to 
model  2-D  cloud  collisions  with  high  resolution.  Self-gravity  was  not  included,  but  a 
range  of  thermodynamic  processes  was.  The  NTSI  arose  when  the  clouds  were  seeded 
with  surface  perturbations  prior  to  collision.  Models  without  surface  perturbations 
did  not  show  the  NTSI,  which  lends  support  to  the  notion  that  the  instability  is 
strictly  nonlinear.  When  radiative  cooling  was  included  in  their  models,  KW  pro- 
duced results  similar  to  the  previous  computational  models.  They  also  considered 
isothermal  models  in  order  to  compare  the  results  with  Vishniac’s  (1994)  analytical 
theory.  They  found  that  a higher  grid  resolution  was  needed  to  see  the  NTSI  in  the 
isothermal  case.  The  isothermal  NTSI  showed  a similar  morphology  to  the  NTSI 
in  models  with  cooling,  but  the  global  quantities  of  the  post-collision  clouds  varied 
between  the  two  cases.  For  example,  the  isothermal  case  overestimated  the  compres- 
sion in  the  collision  mid-plane.  This  indicates  that  the  cooling  mechanism  used  by 
KW  was  rather  weak.  A strongly  cooling  shock  leads  to  a higher  compression  than 
an  isothermal  shock.  Also,  KW  used  a magnetic  pressure  term  to  limit  the  compres- 
sion in  their  radiative  models,  but  not  in  their  isothermal  models.  The  experimental 
growth  rate  found  by  KW  for  their  isothermal  models  agreed  with  the  prediction  by 
Vishniac  (1994).  KW  also  considered  adiabatic  cloud  collisions.  They  found  that  a 
minimum  compression  of  10  was  needed  to  drive  the  NTSI.  An  adiabatic  compression 
of  10  in  the  limit  of  very  high  Mach  number  (velocity  divided  by  sound  speed)  cor- 
responds to  a value  of  7 = 1.22.  This  result  for  adiabatic  cloud  collisions  is  contrary 
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to  the  earlier  studies  which  emphasized  the  need  for  a strong  cooling  mechanism  to 
drive  the  NTSI. 

Miniati  et  al.  (1997)  considered  cloud  collisions  and  explored  the  difference  in 
behavior  between  isothermal  clouds  and  clouds  with  radiative  cooling.  They  found 
that  symmetric  cloud  collisions  in  the  isothermal  limit  behaved  in  a manner  similar 
to  adiabatic  collisions  and  lead  to  the  disruption  of  merged  clouds.  When  radiative 
cooling  was  applied,  symmetric  clouds  coalesced  into  a single  cloud. 

The  behavior  of  a shock  bounded  slab  and  the  development  of  the  NTSI  has 
been  studied  in  some  detail  in  a variety  of  physical  regimes.  Here,  I revisit  the 
methodology  of  HSWK  with  CFDLIB  and  in  light  of  the  studies  which  have  come 
after.  I examine  the  stability  of  a cold,  dense  slab  formed  by  the  collision  of  identical 
supersonic  gas  streams.  I consider  in  detail  the  effects  of  varying  the  cooling  rates 
on  the  growth  of  the  NTSI  and  the  fragmentation  of  the  slab.  At  issue  is  whether 
the  collision  of  gas  streams  and  the  subsequent  instabilities  give  rise  to  collapsing 
clumps. 

4.2  Models 

The  physical  parameters  chosen  for  these  models  are  characteristic  of  cold 
molecular  cloud  cores.  The  gas  is  treated  as  molecular  hydrogen  with  a molecular 
weight  //  = 2 and  a ratio  of  specific  heats  7 = 1.4.  Low  resolution  (129x129  cells) 
models  are  assigned  a cell  size  oi  dx  — dy  = 5 x 10^^  cm.  The  total  spatial  extent  of 
the  grid  is  therefore  6.45  x 10^^  cm,  or  about  0.21  pc.  Observations  of  cloud  cores  with 
sizes  of  the  order  of  0.1  pc  reveal  typical  kinetic  temperatures  of  ~10  iL  and  number 
densities  ~10'*cm“^  (e.g.  Lada  et  al.  1993  and  references  therein).  Accordingly,  an 
initial  density  of  pi  = 10“^^  g cm~^  is  used,  which  corresponds  to  a number  density  of 
n(H2)  ~ 300.  The  initial  temperature  is  Ti  = lOOTL,  but  efficient  cooling  mechanisms 
(see  section  2.3)  drive  the  temperature  lower  in  both  the  pre-shock  and  post-shock 
regions.  The  pre-shock  velocities  of  the  streams  are  opposite  in  direction  and  of  the 
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same  magnitude,  V{  = ±.bkms~^.  The  streams  are  supersonic.  The  Mach  number 
varies  with  the  temperature,  which  in  turn  depends  upon  the  cooling  rate.  In  the 
absence  of  cooling,  the  Mach  number  for  flows  with  temperatures  of  100  it'  is  M = 6.6. 
But  if  the  gas  cools  to  40  K before  the  reaching  the  shock,  a lower  sound  speed  results 
in  a higher  Mach  number,  M = 10.4.  The  changing  dynamical  behavior  of  the  gas 
as  the  amount  of  cooling  is  varied  is  a central  theme  of  this  study.  The  amount  of 
cooling  is  described  by  the  parameter  q,  defined  as  the  fraction  of  the  full  cooling 
rate  used  in  a simulation.  A model  with  q = 1 applies  the  full  cooling  rate;  whereas, 
a model  with  gf  = 0.3  applies  a cooling  rate  that  is  reduced  to  30%  of  its  full  value. 

Each  simulation  begins  with  the  gas  in  the  lower  half  of  the  grid  streaming 
upward  and  the  gas  in  the  upper  half  of  the  grid  streaming  downward,  with  a static 
zone  in  the  center.  The  width  of  the  static  zone  varies,  but  typical  values  are  three 
or  five  cells.  The  bottom  and  top  of  the  grid  (the  horizontal  boundaries)  have  inflow 
boundary  conditions.  Throughout  the  simulation,  material  enters  the  grid  with  the 
same  density,  temperature,  and  velocity  as  the  initial  flows.  Along  the  left  and  right 
sides  (the  vertical  boundaries),  the  gas  is  allowed  to  flow  into  and  out  of  the  grid  freely. 
At  each  computational  step,  the  total  mass  on  the  grid  increases  by  an  amount  equal 
to  the  inflowing  mass  at  the  horizontal  boundaries,  plus  the  inflowing  mass  from  the 
vertical  boundaries,  minus  the  amount  of  mass  which  escapes  the  grid  at  the  vertical 
boundaries.  The  inflow  and  outflow  along  the  vertical  boundaries  can  result  in  the 
generation  of  vortices,  but  typically  do  not  alter  greatly  the  evolution  of  the  slab. 

When  the  post-shock  material  is  allowed  to  cool,  the  colliding  flows  give  rise 
to  a cold,  dense  slab  which  continues  to  accrete  matter.  The  slab  remains  stable  for 
a long  time  (about  1.6  x 10^^  s)  before  numerical  noise  at  the  boundaries  gives  rise  to 
an  instability  which  propagates  into  the  slab.  Perturbations  are  induced  in  the  slab 
well  before  this  time  to  minimize  the  effects  of  the  noise.  The  perturbation  is  applied 
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to  the  y component  of  the  velocity  and  has  the  form 

Vy  = Vy-\-  O-Olt^i  COS  (4.2) 

where  k is  the  wavenumber  of  the  perturbation,  and  yg  the  position  of  the  slab  center. 
The  slab  is  allowed  to  reach  a high  density  through  accretion  prior  to  being  perturbed. 

For  models  requiring  self-gravity,  the  potential  is  found  using  MUDPACK 
in  a manner  very  similar  to  the  method  used  for  the  SGI  models  in  Chapter  3. 
However,  the  process  is  a bit  more  complicated  due  to  the  mass  inflow.  Periodic 
boundary  conditions  are  applied  along  the  vertical  boundaries,  but  are  not  suitable 
at  the  horizontal  boundaries.  The  horizontal  boundaries  require  conditions  that 
change  as  the  mass  on  the  grid  increases.  For  these  models,  the  boundary  values 
for  the  gradient  of  the  potential  are  approximated  by  considering  the  potential  of  an 
infinitely  extended  plane.  During  each  computational  cycle,  an  average  slab  density 
Ps  is  determined  such  that  the  amount  of  mass  in  a uniform  slab  of  density  Ps  matches 
the  mass  of  the  perturbed  (or  unperturbed)  slab.  The  contribution  to  the  potential 
from  the  pre-shock  flows  is  ignored  due  to  the  low  density  in  the  flows  relative  to 
the  cold  slab.  The  potential  gradient  for  an  infinite  plane  of  density  pg  and  thickness 
Ay  is  a constant  in  space,  d^/dy  = 2'KGpgAy.  For  an  unperturbed  dense  slab,  this 
approximation  works  well.  As  instabilities  grow  and  the  slab  begins  to  fragment, 
the  approximation  fails.  Any  results  from  self-gravitating  Cartesian  models  must  be 
viewed  with  this  limitation  in  mind.  Self-gravity  is  not  applied  from  the  beginning  of 
the  simulation,  but  after  the  dense  slab  is  formed.  A sufficient  amount  of  mass  must 
accumulate  before  gravitational  forces  can  compete  with  the  dynamic  forces. 

4.3  Results 

The  dynamic  instability  of  a slab  and  the  role  played  by  the  thermodynamics 
of  the  gas  are  examined  using  models  without  self-gravity.  Figure  (4.1)  shows  the 
growth  of  an  accretion  bounded  slab  formed  by  the  collision  of  adiabatic  gas  streams. 
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The  post-shock  temperature  is  much  higher  than  the  pre-shock  temperature.  Internal 
thermal  pressure  supports  the  slab  against  the  ram  pressure  of  the  accreting  flows. 
The  hot,  dense  slab  quickly  grows  in  width.  The  compression  of  a shock  formed  by 
adiabatic  streams  with  Mach  number  M is  (Harlow  and  Amsden  1971) 

El  = l±i (4  31 

where  Ps  is  the  post-shock  slab  density  and  pi  the  pre-shock  density.  For  flows  of  very 
high  Mach  number,  the  compression  approaches  an  upper  limit  of  (7-I- 1)/(7  — 1).  For 
a molecular  gas  with  7 = 1.4,  the  maximum  compression  is  6.  The  contours  in  figure 
(4.1)  show  that  compression  at  the  center  of  the  slab  is  greater  than  6.  This  is  better 
seen  in  the  line  cuts  shown  in  figure  (4.2).  The  density,  temperature,  pressure,  and  y- 
velocity  are  plotted  along  a constant  value  of  x at  t = 0.8  x lO^^s.  The  grid  resolution 
is  insufficient  to  resolve  the  initial  shock.  As  a result,  the  density  profile  has  a central 
peak,  and  the  temperature  profile  has  a central  dip.  The  two  cancel  each  other  out  in 
determining  the  pressure,  which  remains  flat  inside  the  slab.  Such  behavior  at  shocks 
is  common  among  grid  codes  which  attempt  to  take  numerical  derivatives  across  a 
discontinuity.  The  density  spike  can  be  reduced  by  applying  additional  temperature 
and  momentum  diffusion  to  smear  the  shock  more.  However,  doing  so  results  in  a 
slab  more  prone  to  fragmentation  when  cooling  is  applied.  Over  time,  the  central 
density  peak  decreases  in  the  adiabatic  case.  For  models  with  molecular  cooling, 
the  temperature  inside  the  slab  falls  far  below  the  dip  seen  for  the  adiabatic  case. 
The  central  density  spike  at  the  point  of  collision  serves  only  to  enhance  the  overall 
cooling  by  a small  amount.  Its  presence  has  no  bearing  on  the  evolution  of  cooling 
models. 

When  radiative  cooling  is  added  to  the  simulations,  the  slab  grows  to  a far 
greater  density  than  is  possible  for  an  adiabatic  shock.  The  maximum  slab  density 
attainable  is  constrained  in  each  model  by  the  amount  of  time  that  passes  before 
instabilities  arise  in  the  unperturbed  slab.  Examples  are  shown  in  figures  (4.3)  and 
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Figure  4.1:  Propagation  of  an  adiabatic  7 = 1.4  accretion  shock.  Density  contours 
are  labeled  in  units  of  the  initial  density  (p^  = 10”^^p  cm~^)  and  shown  at  times  of 
a)0.4  X lO^^s;  b)0.8  x lO^^^.  c)1.2  x lO^^^.  d)i.6  x lO^^^.  e)2.0  x lO^^^.  ^^d  f)2.4  x lO^^^ 
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Figure  4.2:  Adiabatic  accretion  shock  line  cuts  at  t = 0.8  x lO^^s.  Cuts  along  a 
constant  value  of  x show  the  a)density;  b) temperature;  c)pressure;  and  d)y- velocity. 


(4.4)  for  models  with  q = 1 (full  cooling)  and  q = 0.3,  respectively.  The  q = I model 
maintains  slab  stability  a short  while  longer,  but  otherwise  the  two  results  are  very 
similar.  In  both  cases,  the  slab  remains  flat  and  stable  until  t ~ 1.6  x lO^^s,  beyond 
which  wiggles  are  seen  near  the  vertical  boundaries.  In  order  to  insure  that  boundary 
effects  do  not  modify  the  growth  of  instabilities  in  the  slab,  an  induced  perturbation 
is  applied  long  before  the  slab  would  otherwise  go  unstable.  Perturbation  times 
of  5 X lO^^s  and  lO^^s  are  considered.  Inflows  into  the  slab  along  the  vertical 
boundaries  result  in  slightly  higher  densities  at  the  edges  than  is  found  throughout 
the  rest  of  the  slab.  The  maximum  slab  density  before  instability  takes  hold  is  about 
3000pi.  This  is  a much  higher  compression  than  that  seen  by  HSWK,  who  noted  a 
maximum  compression  of  about  700pi  for  inflows  of  bkms~^.  They  attributed  their 
low  compression  to  the  high  numerical  diffusion  inherent  to  their  code.  KW  included 
a magnetic  pressure  term  in  their  models  which  limited  the  compression  to  a factor 
of  50.  In  contrast,  the  clond  collision  models  of  Miniati  et  al.  (1997),  which  included 
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Figure  4.3:  Instability  development  for  an  unperturbed  slab  with  q = 1.  Density 
contours  are  labeled  in  units  of  the  initial  density  {pi  = cm“^)  and  shown  at 

times  a)1.5  x lO^^s;  b)1.6  x lO^^s;  c)1.7  x lO^^s;  d)1.8  x lO^^s;  e)1.9  x lO^^s;  and 
f)2.0  X lO^^s.  The  maximum  compression  in  the  center  of  the  slab  is  labeled  Ps-  Note 
how  instabilities  form  at  the  edges  of  the  grid  before  propagating  towards  the  center. 


106 


t=1.7xlO'-’s 

p,=3238pi 


-60  -*0 


X cell 
(c) 


t=1.9x10’-’s 

p,=3630pi 


0 


contour 

0=4567 

b=3262 

c=1958 

d=653. 


I ' ' ' I 


contour 

0=4583 

b=3273 

c=1964 

d=655. 


X can 
(0 


Figure  4.4:  Instability  development  for  an  unpertnrbed  slab  with  q = 0.3.  Density 
contours  are  labeled  in  units  of  the  initial  density  {pi  = gcm~^)  and  shown  at 
times  a)1.5  x lO^^s;  b)1.6  x lO^^s;  c)1.7  x lO^^s;  d)1.8  x lO^^s;  e)1.9  x lO^^s;  and 
f)2.0  X lO^^s.  The  maximum  compression  in  the  center  of  the  slab  is  labeled  p*. 
Instabilities  form  at  the  edges  and  propagate  towards  the  center. 
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Figure  4.5:  Temperature  profiles  for  unperturbed  models  at  1.5  x lO^^s.  Each  plot 
shows  a line  cut  along  constant  x for  a cooling  model  with  a)q  = 1;  and  h)q  = 0.3. 

radiative  cooling,  reached  a compression  of  about  10^.  The  NTSI  shows  similar 
behavior  in  all  of  these  studies,  which  indicates  that  the  strength  of  compression  is 
not  critical  provided  it  exceeds  some  minimum  value.  Indeed,  results  presented  below 
demonstrate  that  the  instability  and  fragmentation  of  the  slab  proceed  in  a nearly 
identical  fashion  whether  the  slab  density  at  the  time  of  perturbation  is  ~ lOOOpi  or 

~ 2000pi. 

The  initial  temperature  of  the  inflowing  gas  is  100  K.  The  flows  begin  to  cool 
once  they  are  within  50  cells  (2.5  x lO^^cm)  of  the  slab.  The  temperature  of  the 
gas  as  it  strikes  an  accretion  shock  depends  upon  the  value  of  q.  Figure  (4.5)  shows 
the  temperature  profiles  at  1.5  x lO^^s  for  the  models  shown  in  figures  (4.3)  and 
(4.4).  When  full  cooling  is  applied  {q  = 1),  the  gas  cools  to  iO  K before  reaching 
the  accretion  shock.  The  temperature  of  the  shocked  gas  rises  to  129  K before  the 
high  slab  density  drives  the  gas  to  cool  to  the  minimum  of  b K.  In  a model  with  a 
lower  cooling  rate  {q  = 0.3),  the  gas  cools  to  12  K before  reaching  the  shock,  spikes 
to  a temperature  of  348  K,  and  then  quickly  cools  to  5 K.  Despite  the  difference  in 
temperature  jump  at  the  shock  between  the  two  models,  the  cooling  is  sufficiently 
strong  in  both  to  lower  the  temperature  to  the  lowest  value  allowed.  As  a result,  the 
central  slab  densities  in  the  models  are  very  similar.  Only  when  perturbations  are 
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imposed  upon  the  slab  are  significant  structural  differences  seen  among  models  with 
varying  degrees  of  cooling. 

4.3.1  Earlier  Perturbation  Time 

When  a velocity  perturbation  (equation  4.2)  is  imposed  upon  the  slab,  mod- 
els with  different  cooling  rates  produce  very  different  results.  The  perturbation  is 
applied  at  t = 5 x lO^^s  when  the  density  in  the  center  of  the  slab  is  ps  = 945p,. 
For  the  models  shown  here,  three  full  waves  are  placed  across  the  grid,  resulting  in  a 
wavenumber  of  k = 2tt/X  = 2.92  x 10~^^cm“^.  Figures  (4.6-4.12)  show  the  instability 
growth  for  perturbed  slabs  in  models  utilizing  different  cooling  rates.  The  base  10 
logarithm  of  the  density  is  plotted  to  highlight  the  morphology  of  the  evolved  slabs. 
The  results  demonstrate  that  the  thermodynamic  properties  of  the  ISM  have  a pro- 
found affect  upon  the  evolution  of  a perturbed  slab.  When  the  full  cooling  rate  is  used 
(figure  4.6)  the  perturbation  develops  into  well  defined  extrema.  Velocity  streaming 
inside  the  slab  drives  the  formation  of  Kelvin- Helmholtz  instabilities  which  give  rise 
to  secondary  peaks.  Secondary  peaks  between  the  extrema  quickly  grow,  so  that  the 
initial  signature  of  the  perturbation  is  largely  erased.  In  contrast,  models  with  much 
lower  cooling  rates,  figures  (4.10-4.11),  develop  extrema  more  slowly  and  preserve 
the  shape  of  the  perturbation  for  longer  times.  Kelvin-Helmholtz  instabilities  devel- 
op more  slowly  or  not  at  all.  There  is  a lower  limit  to  this  trend.  The  ^ = 0.1  model 
shown  in  figure  (4.12)  develops  slowly  in  a way  quite  unlike  all  other  models.  This 
indicates  that  there  exists  a lower  limit  to  the  amount  of  cooling  required  to  drive 
the  NTSI.  All  of  the  models  share  the  characteristic  that  the  density  range  changes 
very  little  once  the  instability  has  grown  to  an  appreciable  size.  The  inflowing  gas 
acts  to  further  displace  and  enlarge  the  slab  rather  than  compress  it.  As  the  pertur- 
bation grows,  the  pre-shock  gas  cools  by  different  amounts  before  reaching  the  shock, 
resulting  in  a temperature  gradient  on  the  shock  surface.  In  most  cases  (the  q = 0.2 
case  being  clearest  exception)  the  instability  growth  is  suppressed  along  the  vertical 


109 


Figure  4.6:  Instability  development  for  a perturbed  slab  with  q = 1.  Logio  density 
contours  are  shown  at  times  a)7.0  x lO^^s;  b)8.0  x lO^^s;  c)9.0  x lO^^s;  d)1.0  x lO^^s; 
e)l.l  X lO^^s;  and  f)1.2  x lO^^s.  The  time  of  perturbation  is  5.0  x lO^^s. 
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Figure  4.7:  Instability  development  for  a perturbed  slab  with  q = 0.8. 
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Figure  4.8:  Instability  development  for  a perturbed  slab  with  q = 0.6. 
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Figure  4.9;  Instability  development  for  a perturbed  slab  with  q = 0.4. 
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Figure  4.10:  Instability  development  for  a perturbed  slab  with  q = 0.3. 
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Figure  4.11:  Instability  development  for  a perturbed  slab  with  q = 0.2. 
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Figure  4.12:  Instability  development  for  a perturbed  slab  with  q = 0.1. 


116 


boundaries  where  inflows  and  outflows  are  allowed.  The  q = 0.2  and  q = 0.3  models 
resemble  the  low  resolution  results  of  KW.  The  models  with  an  intermediate  q value 
appear  to  correspond  to  the  results  of  HSWK,  but  a direct  comparison  is  difficult 
since  their  models  did  not  include  induced  perturbations. 

For  a more  detailed  analysis,  I concentrate  on  two  models  which  highlight  the 
differences  between  a fast  cooling  rate  (the  full  cooling  q = 1 model)  and  a slow 
cooling  rate  (the  q = 0.3  model).  First,  I consider  the  model  with  full  cooling.  In 
figure  (4.13),  the  density  distribution  at  t = 1.1  x lO^^s  is  recreated  as  a grayscale 
image  with  velocity  vectors  superimposed.  (Note  that  the  density  scale  is  opposite 
to  that  used  in  chapter  3.  Here  the  darkest  regions  are  the  densest.)  The  largest 
density  of  p = 1300pi  is  found  along  the  left  boundary.  Excluding  the  boundaries, 
the  maximum  density  in  the  slab  is  505pi,  which  is  roughly  one  half  the  slab  density 
at  the  time  of  perturbation.  The  lowest  density  inside  the  slab  is  0.46pj.  Despite  the 
continuing  inflow,  the  internal  density  falls  as  the  slab  grows  in  size  and  spreads  its 
mass  over  a larger  volume.  For  an  unperturbed  slab,  the  density  after  1.1  x lO^^s  is 
2120pi.  The  instability  results  in  a much  lower  density  than  that  found  in  a stable 
accreting  slab.  The  velocity  vectors  show  that  the  impacting  flows  form  channels 
which  penetrate  through  most  of  the  slab  before  encountering  material  streaming 
from  the  other  side.  Dense  filaments  pervade  the  slab  and  surround  evacuated  regions. 
In  the  most  tenuous  regions  inside  the  slab,  the  gas  is  nearly  static.  Evidence  of 
vortices  is  seen,  but  the  bulk  of  the  flow  maintains  a straight  path  through  the  slab. 

The  q = 0.3  model  evolves  very  differently,  yet  retains  many  of  the  same 
characteristics.  Figure  (4.14)  shows  the  density  distribution  and  velocity  vectors  at 
t = 1.1  X lO^^s.  Once  again,  the  largest  density,  p — 1010 pi,  appears  at  an  edge. 
The  densest  region  away  from  the  boundaries  has  a density  of  A73pi.  Nowhere  inside 
the  slab  does  the  density  fall  below  the  inflow  density.  As  in  the  q — I model,  the 
accreting  flows  strike  a peak  in  the  displacement  and  are  channeled  along  the  surface 
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Figure  4.13:  Perturbed  slab  morphology  for  a g = 1 model  att  = 1.1  x lO^^s.  Velocity 
vectors  are  superimposed  upon  a grayscale  image  of  logio(p)-  Only  vectors  with  a 
magnitude  below  2 x lO^cms"^  are  plotted.  Regions  not  represented  with  vectors 
have  flows  of  higher  magnitude. 
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Figure  4.14:  Perturbed  slab  morphology  for  a g = 0.3  model  at  t = 1.1  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x 10^ cm  are  plotted. 
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and  into  the  slab.  The  channels  for  the  q = 0.3  model  are  more  widely  spaced  due 
to  the  lack  of  secondary  peaks.  Internal  filaments  and  regions  of  vorticity  are  clearly 
seen  inside  the  slab.  When  the  q = 0.3  model  is  perturbed  using  a sine  function  rather 
than  a cosine  function,  the  results  are  similar.  Figure  (4.15)  shows  the  density  and 
velocity  vectors  at  t = 1.1  x lO^^s.  Only  velocity  magnitudes  below  10^cms~^  are 
plotted  in  order  to  bring  out  the  low  velocity  vortices  in  the  slab  interior.  Note  the  well 
defined  dense  and  tenuous  regions  inside  the  slab.  The  velocity  vectors  show  regions 
of  vorticity  adjacent  to  the  dense  filaments,  eis  well  as  oppositely  directed  streams 
bordering  the  tenuous  regions.  Oppositely  directed  streams  may  give  rise  to  Kelvin- 
Helmholtz  instabilities  within  the  slab.  The  lobes  along  the  vertical  boundaries  are 
large  regions  of  circulation. 

The  addition  of  self-gravity  to  the  simulations  does  not  greatly  affect  the 
results,  although  slab  morphologies  are  altered.  For  both  the  q = 1 and  q = 0.3 
models,  self-gravity  was  turned  on  at  the  same  time  the  slab  was  perturbed.  The 
slab  density  at  t = 5 x lO^^s  is  about  945pi,  which  corresponds  to  a Jeans  length 
(equation  2.38)  of  1.20  x lO^^cm  for  a 5 /T  gas.  This  is  smaller  than  the  2.15  x lO^^cm 
perturbation  wavelength.  However,  since  most  of  the  mass  on  the  grid  is  contained 
inside  the  slab,  a better  criterion  for  stability  is  the  minimum  wavelength  for  the 
instability  of  a uniform  sheet.  For  a sheet  with  a surface  density  E,  the  critical 
wavelength  for  gravitational  collapse  is  (Binney  and  Tremaine  1987) 

')RT 


\r  = 


(4.4) 


Using  E = PsAx,  the  result  for  a 5 K gas  is  Ac  = 9.23  x lO^^cm.  This  exceeds 
the  total  spatial  extent  of  the  grid  and  is  more  than  four  times  the  perturbation 
wavelength.  The  slab  is  not  unstable  to  global  collapse.  For  the  q = I model  very 
little  difference  is  seen  between  gravitating  and  non-gravitating  models.  Figure  (4.16) 
shows  the  density  distribution  and  velocity  vectors  at  t = 1.1  x lO^^s  for  the  q — I 
case.  The  peaks  are  somewhat  sharper  and  gravity  has  suppressed  the  number  of 
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Figure  4.15:  Morphology  for  a g = 0.3  model  perturbed  with  a sine  function.  Only 
velocity  vectors  with  a magnitude  below  1 x lO^cm  are  plotted. 
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Figure  4.16:  Morphology  of  a self-gravitating,  perturbed  slab  with  q = 1 aX  t = 
1.1  X lO^^s.  Only  velocity  vectors  with  a magnitude  below  2 x 10^ cm  are  plotted. 
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secondary  peaks  along  the  upper  surface  of  the  slab.  The  largest  density  of  1190pi  is 
lower  than  the  maximum  observed  for  the  non-gravitating  model  and  is  found  along 
the  right  edge  of  the  grid.  The  largest  density  away  from  the  boundaries  is  GOOpi,  a 
value  larger  than  that  for  the  non-gravitating  case.  The  q = 0.3  case  shows  a greater 
difference  between  the  gravitating  and  non-gravitating  cases.  Figure  (4.17)  shows  the 
density  distribution  and  velocity  vectors  at  t = 1.1  x lO^^s  for  the  q = 0.3  case.  The 
shape  of  the  perturbation  is  once  again  preserved,  but  additional  ripples  have  formed 
along  the  slab  surface.  The  flow  pattern  inside  the  slab  has  been  altered  as  well.  The 
maximum  density  along  the  left  edge  is  1740^,.  Away  from  the  edges  the  maximum 
density  is  653pi.  The  q = 0.3  model  produces  regions  of  higher  density  than  the  q = 1 
model  when  self-gravity  is  included.  However,  in  neither  case  are  any  clumps  seen 
of  sufficient  density  to  undergo  collapse.  The  slab  instability  always  results  in  lower 
densities  than  those  in  unperturbed  slabs  allowed  to  continuously  accrete  material. 

4.3.2  Later  Perturbation  Time 

The  same  experiments  were  repeated  for  the  q = 1 and  q = 0.3  models  using 
a later  perturbation  time.  The  form  and  amplitude  of  the  perturbation  are  the  same. 
The  new  time  of  perturbation  is  t = lO^^s,  so  the  slab  is  allowed  to  accrete  material  for 
twice  as  long  as  in  the  previous  models.  The  slab  density  at  the  time  of  perturbation 
is  Ps  = 1920pi.  (As  previously  noted,  along  the  vertical  boundaries  the  density  of 
the  unperturbed  slab  is  higher.)  The  slab  evolution  for  the  q = 1 model  (figure  4.18) 
is  quite  different  than  when  perturbed  earlier.  Internal  flows  are  harder  to  initiate 
inside  the  denser  slab,  so  Kelvin-Helmholtz  instabilities  develop  more  slowly.  Fewer 
secondary  peaks  arise,  and  at  late  times  the  initial  form  of  the  perturbation  can  still 
be  distinguished.  The  density  distribution  is  shown  in  more  detail  along  with  the 
velocity  vectors  in  figure  (4.19)  for  t = 1.6  x lO^^s.  The  highest  density  within  the 
slab  (but  away  from  the  boundaries)  is  1030pj,  which  is  again  about  one-half  the 
density  at  the  time  of  perturbation.  The  minimum  density  within  the  slab  is  0.70pj. 
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Figure  4.17:  Morphology  of  a self-gravitating,  perturbed  slab  with  q = 0.3  at  t 
1.1  X lO^^s.  Only  vectors  with  a magnitude  below  2 x 10^ cm  s ^ are  plotted. 
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Figure  4.18:  Evolution  oi  a.  q = 1 model  perturbed  at  1.0  x lO^^s.  Logio  density 
contours  are  shown  at  times  a)1.2  x lO^^s;  b)1.3  x lO^^s;  c)1.4  x lO^^s;  d)1.5  x lO^^s; 
e)1.6  X lO^^s;  and  f)1.7  x lO^^s. 
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Figure  4.19:  Perturbed  slab  morphology  for  a g = 1 model  at  t = 1.6  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x 10^ cm  are  plotted. 
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Inside  the  slab,  the  densest  patches  of  gas  are  filaments  formed  by  the  collision  of 
oppositely  directed  channels.  This  is  in  contrast  to  the  earlier  perturbation  model 
in  which  the  densest  filaments  are  formed  by  a single  upward  or  downward  directed 
stream. 

The  evolution  of  the  slab  perturbed  at  t = lO^^s  for  the  q = 0.3  case  is  shown 
in  figure  (4.20).  This  model  shows  less  variation  with  the  change  in  perturbation  time 
than  the  q = 1 model.  Overall,  the  morphology  of  the  evolved  slab  remains  largely 
unchanged.  One  difference  is  that  less  structure  is  seen  throughout  the  slab  interior 
as  the  instability  grows  than  in  the  earlier  perturbation  model  (figure  4.10).  Figure 
(4.21)  shows  the  density  distribution  and  velocity  vectors  at  t = 1.6  x lO^^s.  Large 
vortices  appear  in  the  slab  interior.  The  maximum  density  along  the  boundaries  is 
1330/?i,  and  the  densest  region  away  from  the  boundaries  is  not  much  lower  at  1170pi. 
This  density  is  located  on  the  right-most  peak  at  (45,-35).  It  is  interesting  to  note 
that  if  the  cells  along  the  boundaries  are  excluded,  the  q — 1 model  produces  higher 
densities  than  the  q = 0.3  model  when  the  slab  is  perturbed  at  t = 5 x lO^^s,  but  the 
opposite  occurs  when  the  slab  is  perturbed  at  lO^^s. 

The  addition  of  self-gravity  to  the  simulations  is  expected  to  have  a larger 
significance  when  the  slab  is  perturbed  at  a later  time.  As  before,  self-gravity  is 
turned  on  at  t = 5 x lO^^s.  At  t = lO^^s,  the  density  in  the  unperturbed  slab  is 
Ps  = 1920pi.  From  equation  (4.4),  the  maximum  stable  wavelength  is  4.54  x lO^^cm, 
which  is  still  larger  than  the  perturbation  wavelength.  The  slab  is  stable  to  global 
collapse.  Despite  the  denser  slab,  the  self-gravitating  models  are  not  very  different 
from  the  non-gravitating  models.  Figure  (4.22)  shows  the  gravitating  q = 1 case 
at  t = 1.6  X lO^^s.  Other  than  a change  in  the  location  of  some  filaments  and 
some  alterations  to  the  peaks,  the  morphology  matches  that  of  the  non-gravitating 
case  (figure  4.19).  However,  the  densest  filament  in  the  gravitating  case  (located 
at  (50,0))  has  a density  of  1540pj.  This  is  much  denser  than  any  filament  in  the 
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Figure  4.20:  Evolution  of  a g = 0.3  model  perturbed  at  1.0  x lO^^s.  Logio  density 
contours  are  shown  at  times  a)1.2  x lO^^s;  b)1.3  x lO^^s;  c)1.4  x lO^^s;  d)1.5  x lO^^s; 
e)1.6  X lO^^s;  and  f)1.7  x lO^^s. 


Y cell 


128 


-60  -40  -20  0 20  40  60 

X cell 


Figure  4.21:  Perturbed  slab  morphology  for  a 5 = 0.3  model  at  t = 1.6  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x 10^c?ns“^  are  plotted. 
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non-gravitating  case,  but  still  well  below  the  density  needed  to  initiate  collapse. 
A comparison  between  gravitating  (figure  4.23)  and  non-gravitating  (figure  4.21) 
simulations  for  the  q = 0.3  model  reveals  very  similar  morphologies.  The  only  major 
differences  between  the  two  results  are  the  locations  and  magnitudes  of  the  densest 
regions.  In  the  gravitating  case,  the  highest  density  of  1720pi  is  found  along  the  left 
boundary.  The  maximum  density  inside  the  slab  is  916pi,  a value  lower  than  that  for 
the  non-gravitating  case.  Contrary  to  expectations,  for  the  q = 0.3  case  self-gravity 
has  a weaker  effect  upon  the  growth  of  the  instability  when  the  slab  is  perturbed  at 
a later  time.  This  may  be  an  indication  of  the  inaccuracy  of  the  gravity  algorithm 
when  applied  to  inflow  models  in  Cartesian  geometry. 

In  order  to  examine  the  behavior  of  perturbed  adiabatic  slabs,  a series  of  non- 
cooling models  was  run  with  varying  values  of  7.  In  addition  to  7 = 1.4,  values  of  1.1, 
1.2,  and  1.3  were  considered.  In  each  case,  the  center  of  the  hot  slab  was  perturbed 
before  the  thickness  of  the  slab  was  large.  According  to  the  results  of  KW,  the  NTSI 
appears  in  adiabatic  shocks  provided  the  compression  is  greater  than  ten.  By  this 
criterion,  the  7 = 1.1  and  7 = 1.2  models  should  be  unstable  to  the  NTSI,  while  the 
other  two  are  stable.  My  results  do  not  support  this  idea.  No  instability  was  seen  in 
any  of  the  adiabatic  models.  The  material  inside  the  adiabatic  slab  is  nearly  static, 
so  an  induced  perturbation  does  lead  to  the  displacement  of  material.  However,  the 
displacement  does  not  grow  over  time.  The  same  phenomenon  is  seen  in  all  adiabatic 
models,  regardless  of  7.  Either  the  resolution  is  too  low  in  these  models  to  reproduce 
the  findings  of  KW,  or  the  differences  between  the  physical  processes  (colliding  flows 
versus  colliding  clouds)  lead  to  different  behavior.  Continuous  inflow  in  the  colliding 
flow  models  causes  an  adiabatic  slab  to  rapidly  grow  in  thickness.  The  collision  of 
clouds  of  finite  extent  is  different  in  that  the  amount  of  accreting  material  (and  hence 
the  growth  of  an  adiabatic  slab)  is  finite.  It  is  possible  that  a continuous  inflow  of 
matter  wipes  out  the  NTSI  in  adiabatic  models. 
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Figure  4.22:  Morphology  of  a self-gravitating,  perturbed  slab  with  g = 1 at  t = 
1.6  X lO^^s.  Only  velocity  vectors  with  a magnitude  below  2 x lO^cms”^  are  plotted. 
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Figure  4.23:  Morphology  of  a self-gravitating,  perturbed  slab  with  q = 0.3  at  t = 
1.6  X lO^^s.  Only  velocity  vectors  with  a magnitude  below  2 x lO^cm  s ^ are  plotted. 
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4.4  Discussion 

The  main  conclusion  to  be  reached  from  the  Cartesian  models  is  that  the 
cooling  rate  plays  a controlling  role  in  determining  the  evolution  of  a perturbed  cold 
slab.  In  models  with  strong  cooling,  a perturbed  slab  is  prone  to  the  development 
of  short  wavelength  modes  which  override  the  initial  perturbation.  Weakly  cooling 
models  tend  to  preserve  the  perturbation  wavelength  for  long  times.  In  all  cases,  the 
slab  density  decreases  as  the  gas  is  spread  over  a larger  volume.  Streaming  motions 
within  the  slab  give  rise  to  filaments  and  vortices.  The  slab  ultimately  fragments 
into  a few  dense  islands  of  material,  but  no  collapsing  clumps  are  observed.  Altering 
the  cooling  rate  leads  to  a larger  change  than  adding  self-gravity  to  the  simulations. 
This  may  not  hold  true  for  models  which  have  either  a higher  density  infiow  or  a 
longer  perturbation  wavelength  which  enhances  self-gravity.  Additional  models  with 
different  initial  density  values  would  be  informative.  In  their  cylindrical  models, 
HSWK  found  that  very  massive  slabs  evolved  in  a manner  similar  to  less  massive 
slabs.  The  current  results  support  this  claim  and  can  be  further  tested  and  verified 
for  different  cooling  rates.  A study  of  models  with  different  inflow  velocities  would 
also  be  constructive.  With  one  exception,  the  Cartesian  models  presented  in  this 
chapter  use  the  same  perturbation  function.  Simulations  in  which  a sine  function 
replaced  the  cosine  function  in  the  perturbation  equation  showed  the  same  qualitative 
behavior,  but  no  in-depth  analysis  was  performed.  Models  with  strong  cooling  lose 
the  signature  of  the  perturbation  as  the  instability  grows,  so  altering  the  perturbation 
function  should  have  little  effect,  provided  the  wave  number  remains  unchanged.  Two 
dimensional  models  must  ultimately  give  way  to  three  dimensional  models  in  order  to 
accurately  represent  the  systems  under  consideration.  As  work  proceeds,  additional 
physics  must  be  added  as  well,  including  magnetic  fields,  heating  processes,  and 
obliquely  colliding  flows. 
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The  strong  dependence  of  the  slab  instability  on  the  cooling  rate  casts  doubt 
upon  the  utility  of  Vishniac’s  (1994)  analytic  description  of  the  NTSI  in  an  isothermal 
gas.  However,  it  should  be  noted  that  despite  the  many  differences  among  models 
with  varying  cooling  rates,  the  growth  rate  for  the  slab  displacement  is  similar  in  most 
cases.  Attempts  to  find  evidence  for  a linear  growth  phase  to  the  NTSI  have  failed, 
upholding  that  aspect  of  Vishniac’s  theory.  Furthermore,  a higher  wave  number  leads 
to  a faster  growth  rate,  as  expected  from  theory.  Many  authors  have  noted  that  the 
isothermal  limit  is  not  a good  representation  of  interstellar  gas  (e.g.  HSWK,  Vazquez- 
Semadeni  et  al.  1996).  For  a gas  that  can  be  described  as  polytropic,  such  that 
P oc  p'’'®,  the  effective  value  of  7 for  strongly  cooling  flows  is  less  than  unity,  7e  < 1. 
New  theories  are  needed  which  take  into  account  the  high  compressions  afforded  by 
realistic  cooling  rates  in  the  ISM.  The  nonlinear  nature  of  the  instability  makes  the 
development  of  any  such  theory  very  difficult.  Furthermore,  because  different  phases 
of  the  ISM  exhibit  different  heating  and  cooling  processes,  the  behavior  of  the  NTSI 
is  expected  to  vary  in  different  regions  of  the  ISM. 

The  Cartesian  simulations  presented  in  this  chapter  show  many  similarities  to 
previous  simulations  performed  in  cylindrical  geometry.  Models  with  low  cooling  rates 
exhibit  morphologies  very  similar  to  the  models  of  MK,  even  though  the  compression 
in  the  slab  is  much  higher.  I find  no  evidence  of  the  NTSI  in  perturbed  adiabatic 
slabs.  This  may  be  due  to  a lack  of  sufficient  resolution  or  the  result  of  having 
continuously  accreting  shocks.  Further  comparisons  with  previous  work  are  reserved 
for  the  next  chapter  in  which  colliding  flows  in  cylindrical  geometry  are  considered. 


CHAPTER  5 

COLLIDING  FLOWS  IN  CYLINDRICAL  GEOMETRY 

For  studies  of  self-gravitating  systems  using  two  dimensional  modeling,  cylin- 
drical geometry  is  preferable  to  Cartesian  geometry.  Azimuthal  symmetry  imparts 
a three  dimensional  nature  to  the  simulations;  clumps  which  form  along  the  sym- 
metry axis  approximate  spheres.  Previous  computational  simulations  of  the  NTSI 
and  slab  fragmentation  (e.g.  HWSK,  SBP,  KW)  focussed  on  cylindrical  geometry.  In 
order  to  better  compare  new  results  with  the  previous  studies,  cylindrical  models  are 
performed  with  and  without  self-gravity.  The  formation  of  gravitationally  collaps- 
ing clumps  is  investigated.  Five  groups  of  models  are  considered:  models  with  free 
outer  boundaries  and  full  cooling,  models  with  reflective  outer  boundaries  and  full 
cooling,  models  with  free  outer  boundaries  and  reduced  cooling,  models  with  varying 
minimum  temperatures,  and  models  with  a higher  inflow  velocity. 

5.1  Models 

In  order  to  differentiate  the  cylindrical  models  from  the  Cartesian  models,  I 
use  “disk”  rather  than  “slab”  to  describe  the  shock-bounded  object.  The  formation 
and  subsequent  instability  of  disks  and  slabs  are  very  similar.  Colliding  gas  flows 
give  rise  to  a dense  disk  bounded  by  accretion  shocks.  If  allowed  to  cool,  the  disk 
becomes  very  dense  prior  to  the  onset  of  instability.  The  physical  characteristics  for 
the  cylindrical  models  are  similar  to  those  for  the  Cartesian  models.  Unless  otherwise 
noted,  the  inflowing  gas  has  the  values  pt  = cm“^,  Tj  = 100  AT  (prior  to  the 

onset  of  cooling),  and  Vi  = ±5kms~^.  The  cell  size  is  the  same  as  that  for  the 
Cartesian  models,  Ar  = Az  = 5 x lO^^cm.  The  mass  inside  a cylinder  with  a 
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radius  of  97  cells  and  a height  of  129  cells  is  4.77  x 10^^  5'  — 0.24  Mq.  Once  the 
simulation  begins,  mass  enters  the  grid  through  the  horizontal  boundaries  at  a rate 
of  7.39  X 10^°  55“^.  The  symmetry  a.xis  is  represented  by  a reflecting  boundary  on 
the  left  edge  of  the  grid.  The  right  edge  of  the  grid  has  a free  boundary  allowing 
inflows  and  outflows. 

The  application  of  self-gravity  in  cylindrical  coordinates  is  accomplished  using 
the  method  described  in  chapter  2.  Because  of  the  geometry,  gravity  pulls  material 
towards  the  symmetry  axis.  Along  the  symmetry  axis,  the  radial  component  of  the 
potential  gradient  is  zero.  The  boundary  values  of  the  potential  away  from  the  sym- 
metry axis  are  determined  for  each  time  step  by  explicitly  summing  the  contributions 
from  the  annulus  described  by  each  cell.  This  procedure  is  more  accurate  than  the 
method  used  for  self-gravitating  Cartesian  models,  but  also  requires  far  more  com- 
puting time.  The  grid  size  for  most  models  is  kept  rather  small  (97x129  cells)  in 
order  to  keep  the  run  times  reasonable.  Still,  the  grid  has  far  higher  resolution  than 
the  25x50  grid  used  by  HSWK. 

Results  are  presented  for  models  with  two  different  cooling  rates,  g = 1 and 
q = 0.3.  Also,  a number  of  additional  parameter  variations  are  considered  for  the 
q = 1 models.  More  emphasis  is  placed  on  the  q = 1 models  since  they  implement 
the  best  cooling  model  available  for  simulating  molecular  cloud  cores.  As  in  the 
Cartesian  models,  the  cooling  is  confined  to  within  50  cells  of  the  disk.  Simulations 
were  performed  for  a wide  range  of  perturbations,  but  only  a few  are  chosen  for 
detailed  presentation.  In  order  to  examine  the  differences  between  a Cartesian  and  a 
cylindrical  mesh,  simulations  which  have  the  same  perturbation  function  as  the  one 
used  in  chapter  4 (equation  4.2)  are  used.  Other  models  considered  in  this  chapter 
have  a perturbation  of  the  form 


Vz  = Vz  + O.OlUi  sin  (kr)  sin  {nr/rd)e 


(5.1) 


136 


where  is  the  total  radius  of  the  disk  (the  spatial  extent  of  the  grid  in  the  radial 
direction).  The  second  sine  term  superimposes  one-half  of  a sine  wave  (0  < r/r^  < 1) 
upon  the  perturbation  in  order  to  localize  the  perturbation  to  the  center  of  the 
grid  and  minimize  the  initial  growth  along  the  boundaries.  Since  periodic  boundary 
conditions  are  not  used  in  the  gravity  solver,  the  perturbation  function  need  not  be 
strictly  periodic.  A half  integer  number  of  waves  is  excited  across  the  grid  so  that 
both  sine  functions  have  their  maxima  at  the  grid  center. 

5.2  Results 

When  gas  streams  collide  and  are  allowed  to  cool,  a cold,  dense  disk  forms 
that  remains  stable  for  a long  time  if  left  unperturbed.  Figures  (5.1)  and  (5.2)  show 
the  results  of  simulations  performed  on  a 97x  129  grid  for  the  q = \ and  q = 0.3  cases, 
respectively,  for  simulations  performed  on  a 97x129  grid.  (In  all  plots  of  cylindrical 
models,  the  r axis  is  shifted  by  one  cell  so  that  the  symmetry  axis  does  not  lie  along 
the  plot  axis.)  The  disks  remain  stable  until  t ~ 1.4  x lO^^s.  They  do  not  last 
as  long  as  the  Cartesian  models  presented  in  the  previous  chapter.  The  disks  are 
densest  near  the  outer  radial  boundary  due  to  the  inflow  of  material,  but  they  first 
go  unstable  along  the  reflective  symmetry  axis.  Overall,  the  disk  in  the  q = \ model 
is  denser  than  the  disk  in  the  q = 0.3  model.  A stronger  cooling  mechanism  inside 
the  disk  results  in  a lower  thermal  pressure,  allowing  a greater  compression.  A close 
look  at  figures  (5.1)  and  (5.2)  reveals  that  the  onset  of  instability  is  signaled  by  the 
presence  of  high  frequency  ripples  in  the  disks.  For  the  q = I case,  the  ripples  are 
largest  near  the  symmetry  axis  and  tend  to  diminish  with  increasing  r.  The  ripples  in 
the  q = 0.3  case  are  more  uniform  in  size  throughout  the  disk.  Ripples  at  large  r grow 
quickly  enough  to  drive  the  disk  unstable  before  the  instability  from  the  symmetry 
axis  propagates  outward. 

The  switch  from  Cartesian  to  cylindrical  geometry  causes  the  disk  to  evolve 
differently,  even  if  all  else  is  kept  the  same.  The  q — 1 and  q = 0.3  models  from 
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Figure  5.1:  Onset  of  instability  in  an  unperturbed  g = 1 disk.  Logxo  density  contours 
are  shown  at  times  a)1.4  x lO^^s;  b)1.5  x lO^^s;  c)1.6  x lO^^s;  and  d)1.7  x lO^^s.  An 
instability  initially  appears  along  the  symmetry  axis  and  then  propagates  outward. 
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Figure  5.2:  Onset  of  instability  in  an  unperturbed  q = 0.3  disk.  Logio  density 
contours  are  shown  at  times  a)1.4x  lO^^s;  b)1.5x  lO^^s;  c)1.6x  lO^^s;  and  d)1.7 x lO^^s. 
The  disk  first  goes  unstable  along  the  symmetry  axis. 
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the  previous  chapter  are  reproduced  on  129x129  cylindrical  grids  to  demonstrate  the 
changes.  Figure  (5.3)  shows  the  development  of  the  disk  instability  for  a g = 1 disk 
perturbed  at  5.0  x lO^^s.  The  perturbation  function  follows  the  cosine  form  (equation 
4.2)  used  for  the  q = 1 Cartesian  model  (figure  4.6).  The  cylindrical  and  Cartesian 
models  are  similar  in  that  both  show  the  growth  of  secondary  peaks  throughout  the 
disk.  In  the  cylindrical  version,  the  reflective  symmetry  axis  behaves  far  differently 
than  the  outer  boundary.  The  fastest  growing  peak  is  along  the  symmetry  axis.  A 
comparison  of  the  density  values  once  the  disk  is  well  evolved  reveals  additional  dif- 
ferences. The  density  distribution  at  t = 1.1  x lO^^s  is  reproduced  as  a grayscale  plot 
with  overlaying  velocity  vectors  in  figure  (5.4).  Additional  line  contours  are  added 
to  highlight  the  densest  regions.  The  evolved  disk  has  dense  filaments  throughout 
its  interior.  Also  evident  are  several  regions  within  the  disk  in  which  the  density  is 
lower  than  the  initial  density.  The  most  striking  difference  between  the  cylindrical 
and  Cartesian  models  is  the  appearance  of  a dense  filament  along  the  symmetry  axis. 
The  highest  density  in  the  filament  (which  is  also  the  highest  density  anywhere  on  the 
grid)  is  2270 Pi.  This  is  the  first  evidence  of  a density  higher  than  the  density  of  an 
unperturbed  disk  at  the  same  time.  The  density  of  an  unperturbed  slab  at  1.1  x lO^^s 
is  2120pi.  In  the  Cartesian  model  (figure  4.13),  the  largest  observed  density  at  the 
same  time  is  1300pi.  Clearly,  a cylindrical  grid  leads  to  the  geometric  focusing  of 
material  towards  the  symmetry  axis.  Self-gravity  is  expected  to  further  enhance  the 
density  growth  along  the  axis. 

A similar  comparison  is  made  between  cylindrical  and  Cartesian  models  for 
the  q = 0.3  case.  The  results  for  the  cylindrical  model  are  shown  in  figure  (5.5). 
The  cylindrical  and  Cartesian  (figure  4.10)  models  evolve  in  similar  manners.  The 
perturbation  wavelength  is  well  preserved  to  late  times.  As  in  the  cylindrical  q = 1 
case,  the  instability  grows  fastest  along  the  symmetry  axis  in  the  cylindrical  q = 0.3 
model.  However,  the  maximum  density  (p  = 1950pi)  is  found  along  the  outer  radial 
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Figure  5.3:  Instability  growth  for  a g = 1 disk  on  a 129x129  grid.  Logio  density 
contours  are  shown  at  times  a)7.0  x lO^^s;  b)8.0  x lO^^s;  c)9.0  x lO^^s;  d)1.0  x lO^^s; 
e)l.l  X lO^^s;  and  f)1.2  x lO^^s.  The  time  of  perturbation  is  5.0  x lO^^s. 
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Figure  5.4:  Perturbed  disk  morphology  for  a g = 1 model  at  t = 1.1  x lO^^s.  Velocity 
vectors  are  superimposed  upon  a grayscale  image  of  logio(p)-  vectors  with  a 

magnitude  below  2 x lO^cms"^  are  plotted.  Regions  not  represented  with  vectors 
have  flows  of  higher  magnitude.  Line  contours  highlight  the  densest  regions. 
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Figure  5.5:  Instability  growth  for  a g = 0.3  disk  on  a 129x129  grid.  Logio  density 
contours  are  shown  at  times  a)7.0  x lO^^s;  b)8.0  x lO^^s;  c)9.0  x lO^^s;  d)1.0  x lO^^s; 
e)l.l  X lO^^s;  and  f)1.2  x lO^^s.  The  time  of  perturbation  is  5.0  x lO^^s. 
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boundary  rather  than  in  the  filament  which  forms  along  the  symmetry  axis.  Figure 
(5.6)  shows  the  density  distribution  and  velocity  vectors  at  t = 1.1  x lO^^s.  The  long 
filament  along  the  symmetry  axis  has  a maximum  density  of  594pi,  far  lower  than 
that  seen  in  the  q = 1 model  at  the  same  time.  The  maximum  density  is  located  in 
the  same  place  as  seen  in  the  Cartesian  model,  but  the  value  is  much  higher.  The 
tendency  for  material  to  move  towards  the  symmetry  axis  may  be  preventing  some 
material  from  fiowing  out  of  the  grid,  resulting  in  a larger  density  concentration.  As 
in  the  Cartesian  case,  vortices  are  clearly  seen  within  the  interior  of  the  disk. 

5.2.1  Free  Outer  Boundaries  and  Full  Cooling 

Due  to  the  long  run  times  needed  for  the  gravity  solver,  the  bulk  of  the  cylindri- 
cal simulations  were  performed  on  a smaller  97x129  grid.  The  perturbation  function 
contains  two  sine  terms  as  described  by  equation  (5.1).  The  first  set  of  models  has 
a perturbation  time  of  5 x lO^^s.  At  the  time  of  perturbation,  the  density  along 
the  symmetry  axis  is  945pi,  but  increases  at  large  r (figure  5.7)  due  to  inflows  along 
the  outer  boundary.  The  temperature  distribution  follows  the  opposite  trend,  since  a 
higher  density  results  in  more  cooling  and  hence  lower  temperatures.  (The  additional 
inflows  arise  very  early  when  free  boundaries  allow  outflows  and  inflows,  regardless 
of  the  total  extent  of  the  grid.  In  the  next  subsection,  I describe  a series  of  models 
with  reflective  outer  boundaries.  The  overall  results  are  very  similar.)  The  pertur- 
bation wavelength  is  1.39  x lO^^cm;  3.5  full  waves  fit  on  the  grid.  The  results  for  a 
non-gravitating  q = I model  are  shown  in  figure  (5.8).  At  early  times,  the  largest 
growth  is  seen  at  the  center  of  the  grid  where  the  maximum  perturbation  is  applied. 
Secondary  peaks  develop  early,  but  they  do  not  mask  the  perturbation  signature. 
Inflowing  material  is  transferred  along  the  peaks  and  into  the  interior,  giving  rise  to 
dense  filaments  within  the  perturbed  disk.  Some  of  the  material  strikes  a peak  and 
moves  towards  the  vertical  boundaries,  leading  to  outflows  along  the  outer  bound- 
ary. Figure  (5.9)  shows  the  density  distribution  and  velocity  flow  within  the  disk 
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Figure  5.6:  Perturbed  disk  morphology  for  a 5 = 0.3  model  at  t = 1.1  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x 10^cms“^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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Figure  5.7:  Radial  density  and  temperature  distributions  for  an  unperturbed  g = 1 
disk  at  5 X lO^^s. 

at  t = 1.1  X lO^^s.  The  highest  density  (p  = 2930pj)  is  found  in  a filament  along 
the  symmetry  axis.  As  seen  in  previous  q — I models,  the  lowest  density  regions 
are  found  inside  the  disk.  The  internal  fiows  are  largely  uni-directional,  but  some 
vorticity  is  present. 

The  perturbation  wavelength  (A  = 1.39  x lO^^cm)  is  well  below  the  critical 
wavelength  for  the  collapse  of  a uniform  sheet.  From  equation  (4.4),  the  critical  wave- 
length at  the  time  of  perturbation  is  Ac  = 9.23  x lO^^cm  for  p = 945pi  and  T = 5 K. 
A few  additional  critical  values  are  worth  mentioning.  The  total  diameter  of  the  disk 
(accounting  for  symmetry)  is  194  x 5 x lO^^cm  = 9.7  x lO^^cm.  The  minimum  density 
needed  to  drive  the  entire  disk  towards  collapse  if  the  perturbation  wavelength  were 
equal  to  the  disk  diameter  is  pc  = 899 pi.  Clumps  which  form  along  the  symmetry 
axis  can  be  approximated  as  spheres.  From  the  spherical  Jeans  equation  (equation 
2.38),  the  minimum  density  needed  for  catastrophic  collapse  within  a single  cell  is 
Pc  = 5.48  X 19^ Pi-  Collapse  is  assured  should  the  density  reach  this  value  in  any 
cell  along  the  symmetry  axis.  (This  is  a conservatively  high  threshold  for  collapse. 
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Figure  5.8:  Instability  growth  for  a g = 1 disk  perturbed  at  5 x lO^^s.  Logio  density 
contours  are  shown  at  times  a)7.0x  lO^^s;  b)9.0x  lO^^s;  c)l.l  x lO^^s;  and  d)1.3xl0^^s. 
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Figure  5.9;  Morphology  of  a g = 1 disk  at  1.1  x lO^^s.  Only  velocity  vectors  with 
a magnitude  below  2 x 10^cms“^  are  plotted.  Line  contours  highlight  the  densest 
regions. 
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Lower  densities  distributed  across  multiple  cells  can  also  lead  to  collapse.)  For  sim- 
ulations which  include  self-gravity,  the  gravity  solver  is  activated  at  the  time  of  the 
perturbation.  The  density  of  the  disk  at  that  time  is  low  enough  that  gravity  does 
not  have  a substantial  effect  upon  the  instability  right  away.  However,  the  inclu- 
sion of  the  gravity  solver  into  the  code  introduces  an  additional  source  of  numerical 
noise.  The  result  is  an  additional  disturbance  which  begins  at  the  outer  boundary 
and  propagates  towards  the  symmetry  axis.  Figure  (5.10)  shows  the  evolution  of 
this  disturbance.  As  the  gravity  induced  disturbance  propagates  towards  the  axis, 
it  gives  rise  to  a density  spike  on  either  side  of  the  expanded  region  (as  shown  for 
t = 1.1  X lO^^s  in  figure  5.11).  The  disturbance  does  not  reach  the  symmetry  axis 
prior  to  the  onset  of  instability  along  the  axis  observed  in  the  non-gravitating  models. 
The  growth  of  clumps  along  the  axis  is  not  greatly  affected.  The  same  phenomena 
are  observed  for  the  q = 0.3  gravitating  models. 

The  disturbance  which  propagates  inward  in  the  gravitating  models  is  over- 
whelmed when  a perturbation  is  induced  throughout.  The  instability  growth  for 
a self-gravitating  q — I disk  is  shown  in  figure  (5.12).  The  morphology  is  rather 
similar  to  that  of  the  non-gravitating  model  at  early  times,  but  higher  densities  are 
soon  observed  along  the  symmetry  axis.  The  disturbance  seen  in  the  unperturbed 
model  (figure  5.10)  is  still  present,  but  far  denser  clumps  and  filaments  form  along 
the  axis.  Figure  (5.13)  shows  the  density  distribution  and  velocity  flow  within  the 
disk  at  t = 1.1  X lO^^s.  At  this  time,  the  density  along  the  symmetry  axis  has  reached 
a maximum  of  1.43  x 10^ Little  vorticity  is  seen  within  the  disk. 

Unlike  in  the  Cartesian  models,  the  addition  of  self-gravity  to  the  cylindrical 
models  has  a profound  effect.  Without  gravity,  the  perturbed  disk  grows  in  volume, 
and  geometric  focusing  results  in  the  formation  of  dense  filaments  along  the  symmetry 
axis.  The  densities  within  these  filaments  fluctuate  as  clumps  contract  and  re-expand 
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Figure  5.10:  Evolution  of  an  unperturbed  9 = 1 disk  when  gravity  is  activated  at 
5 X lO^^s.  Logio  density  contours  are  shown  at  times  a)7.0  x lO^^s;  b)9.0  x lO^^s; 
c)l.l  X lO^^s;  and  d)1.3  x lO^^s. 


150 


R cell 

(o) 


0 20  40  60  80  100 

R cell 

(b) 


Figure  5.11:  Radial  density  and  temperature  distributions  for  an  unperturbed,  self- 
gravitating  q = 1 disk  at  1.1  x lO^^s. 


in  response  to  the  dynamic  forces  and  their  own  internal  pressure.  Models  with  self- 
gravity are  different  in  that  long-lived,  high  density  concentrations  eventually  form 
along  the  symmetry  axis.  Figure  (5.14)  contains  plots  of  the  maximum  density  along 
the  symmetry  axis  versus  time  for  the  q = \ models.  The  models  represented  in  the 
figure  are  the  unperturbed,  non-gravitating  case  (solid  line),  the  unperturbed,  self- 
gravitating  case  (dotted  line),  the  perturbed,  non-gravitating  case  (dashed  line),  and 
the  perturbed,  self-gravitating  case  (dashed-dotted  line).  The  quantity  logio(p/pi)  is 
plotted  in  order  to  show  the  non-gravitating  and  gravitating  cases  on  the  same  graph. 
A value  of  log^o(p/Pi)  = 5. 74  is  needed  to  reach  the  critical  density  for  the  collapse  of 
a single  cell.  The  unperturbed,  non-gravitating  model  accretes  material  at  a nearly 
constant  rate  until  the  boundary  induced  instability  starts  to  grow  at  t 1.4  X lO^^s. 
The  maximum  density  in  the  perturbed,  non-gravitating  model  oscillates  on  a short 
time  scale,  but  never  deviates  from  the  unperturbed  case  by  more  than  an  order  of 
magnitude.  No  collapsing  clumps  are  observed.  The  gravitating  models  do  produce 
collapsing  clumps  along  the  symmetry  axis,  but  only  after  the  disturbance  induced 
by  the  gravity  solver  reaches  the  axis.  This  happens  earlier  in  the  model  lacking  an 
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Figure  5.12:  Instability  growth  for  a self-gravitating  q = \ disk  perturbed  at  5 x lO^^s. 
Logio  density  contours  are  shown  at  times  a)7.0  x lO^^s;  b)9.0  x lO^^s;  c)l.l  x lO^^s; 
and  d)1.3  x lO^^s.  The  gravity  solver  is  switched  on  at  the  time  of  the  perturbation. 
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Figure  5.13:  Morphology  of  a self-gravitating  q = \ disk  at  1.1  x lO^^s.  Only  velocity 
vectors  with  a magnitude  below  2 x lO^cms"^  are  plotted.  Line  contours  highlight 
the  densest  regions. 
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Figure  5.14:  Maximum  density  along  the  symmetry  axis  versus  time  for  models 
with  a free  outer  boundary.  Represented  are  the  unperturbed,  non-gravitating  case 
(solid  line),  the  unperturbed,  self-gravitating  case  (dotted  line),  the  perturbed,  non- 
gravitating case  (dashed  line),  and  the  perturbed,  self-gravitating  case  (dashed-dotted 
line). 
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induced  perturbation,  as  seen  by  the  sudden  density  rise  at  t ~ 1.6  X lO^^s.  The 
perturbed,  gravitating  case  shows  oscillations  similar  to  those  in  the  perturbed,  non- 
gravitating case  before  reaching  the  point  of  collapse.  An  upward  trend  in  density 
is  observed  starting  at  t ~ 1.4  x lO^^s,  which  would  likely  produce  regions  of  local 
collapse  even  without  the  disturbance  propagating  inward  from  the  outer  boundary. 

The  use  of  cylindrical  geometry  permits  a clear  definition  for  the  total  mass  in 
the  simulation.  For  this  purpose,  the  computational  region  is  taken  to  be  a cylinder 
with  a radius  of  97  cells  {R  — 4.85  x 10^^  cm)  and  a height  of  129  cells  {Z  = 6.45  x 
10^^ cm).  As  matter  enters  the  grid  through  the  top  and  bottom  boundaries,  the 
total  mass  increases  at  a constant  rate  given  by 


m{t)  — TxR^Zpi 


(5.2) 


The  total  mass  will  also  vary  due  to  flows  across  the  outer  radial  boundary.  To 
determine  the  total  mass  at  a given  time,  each  cell  is  treated  as  an  annulus,  in  a 
manner  similar  to  that  used  in  the  gravity  solver.  The  masses  of  all  annuli  are  summed 
to  arrive  at  the  total.  Figure  (5.15a)  shows  the  result  for  the  non-gravitating  9 = 1 
models.  The  change  due  to  inflow  from  the  horizontal  boundaries  is  represented  by 
the  solid  line.  The  total  mass  in  the  unperturbed  model  is  always  higher  than  the 
constant  inflow  value,  indicating  the  presence  of  additional  inflows  along  the  outer 
boundary.  When  a perturbation  is  applied,  some  of  the  incoming  flows  are  deflected 
towards  the  outer  boundary,  resulting  in  a net  outflow.  The  total  masses  for  the 
gravitating  models  as  a function  of  time  are  plotted  in  figure  (5.15b).  As  expected, 
gravity  acts  strongly  to  pull  additional  matter  onto  the  grid. 

5.2.2  Reflective  Outer  Boundaries  and  Full  Cooling 


When  material  is  allowed  to  flow  across  the  onter  bonndary,  an  additional 
inflow  is  established  early  in  the  simulation.  This  additional  inflow  is  responsible  for 
the  rise  in  density  at  large  r shown  in  figure  (5.7),  as  well  as  the  higher  than  expected 
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Figure  5.15:  Total  mass  versus  time  for  models  with  Vi  = bkms~^.  The  total  mass 
is  compared  to  the  mass  due  to  the  constant  inflow  from  the  horizontal  boundaries 
for  a)non-gravitating  models;  and  b)gravitating  models. 
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mass  evident  in  an  unperturbed  disk.  The  disturbance  seen  in  the  gravitating  mod- 
els may  be  the  same  phenomenon  amplified  by  the  effects  of  self-gravity.  In  order 
to  examine  the  effects  of  the  flows  along  the  outer  boundary,  the  simulations  were 
repeated  with  a reflective  boundary  condition  along  the  outer  boundary.  This  creates 
a more  artificial  situation  in  which  no  material  is  allowed  to  flow  across  the  outer 
boundary.  When  left  unperturbed,  the  disk  remains  stable  for  a long  time.  As  in  the 
outflow  case,  the  disk  first  goes  unstable  along  the  symmetry  axis  at  t ~ 1.4  x lO^^s 
(figure  5.16).  The  lack  of  inflows  along  the  outer  boundary  allows  instabilities  to 
develop  there  sooner.  At  the  time  of  perturbation  (t  = 5 x lO^^s),  the  density  and 
temperature  are  constant  throughout  the  entire  disk  (figure  5.17). 

Results  for  a perturbed  disk  are  shown  in  figure  (5.18).  The  instability  growth 
proceeds  in  a very  similar  manner  as  that  in  the  free  boundary  model.  Even  along 
the  outer  boundary  the  two  models  show  very  few  differences.  Away  from  the  outer 
boundary,  the  models  are  nearly  identical.  The  density  distribution  and  velocity  flow 
within  the  disk  are  shown  for  t = 1.1  x lO^^s  in  figure  (5.19).  The  densest  region 
{p  = 2940pi)  is  found  along  the  symmetry  axis.  Switching  from  a free  boundary  to  a 
reflective  boundary  has  very  little  effect  upon  non-gravitating  models. 

Gravitating  models  with  a reflective  outer  boundary  deviate  from  the  gravi- 
tating, free  boundary  models.  The  large  mass  increase  shown  in  figure  (5.15)  for  the 
gravitating  outflow  models  is  prevented  by  the  reflective  outer  boundary.  All  models 
with  a reflective  outer  boundary  have  total  masses  which  increase  at  a constant  rate 
due  to  the  inflows  from  the  horizontal  boundaries.  The  total  mass  never  decreases. 
As  shown  in  figure  (5.20),  an  unperturbed  disk  goes  unstable  sooner  when  gravity 
is  activated.  Instabilities  first  develop  along  both  the  symmetry  axis  and  the  outer 
boundary.  Gravity  pulls  material  towards  the  symmetry  axis,  resulting  in  a growing 
tenuous  region  near  the  outer  boundary. 
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Figure  5.16:  Onset  of  instability  in  an  unperturbed  q = I disk  with  reflective  outer 
boundary.  Logic  density  contours  are  shown  at  times  a)1.4  x lO^^s;  b)1.5  x lO^^s; 
c)1.6  X lO^^s;  and  d)1.7  x lO^^s.  The  disk  first  goes  unstable  along  the  symmetry 


axis. 
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Figure  5.17:  Radial  density  and  temperature  distributions  for  an  unperturbed  q = I 
disk  with  a reflective  outer  boundary  at  5 x lO^^s. 

Results  for  a perturbed,  gravitating  model  are  shown  in  figure  (5.21).  The 
reflective  boundary  model  closely  resembles  the  free  boundary  model  during  the  early 
stages  of  the  instability  growth.  However,  the  two  models  diverge  at  later  times 
when  the  mass  difference  between  them  becomes  considerable.  By  t = 1.1  x lO^^s, 
the  disk  morphology  for  the  reflective  boundary  model  (figure  5.22)  is  substantially 
different  from  that  in  the  free  boundary  model.  The  largest  density  observed  for  the 
reflective  boundary  model  (p  = 3510pi)  is  considerably  less  than  the  largest  density 
observed  for  the  free  boundary  model  (p  = 14,300pj).  However,  the  formation  of 
collapsing  clumps  does  occur  in  the  perturbed,  gravitating  model  with  a reflective 
outer  boundary.  As  shown  in  figure  (5.23),  the  reflecting  boundary  model  reaches 
the  critical  density  earlier  than  the  free  boundary  model.  The  outflows  induced  in 
the  perturbed  free  boundary  model  are  prohibited  in  the  reflective  boundary  model, 
so  an  earlier  collapse  is  possible. 

5.2.3  Free  Outer  Boundaries  and  Reduced  Cooling 

Simulations  were  performed  for  q = 0.3  models  with  a free  outer  boundary. 
Perturbed  models  again  have  a perturbation  time  of  5 x which  is  the  same  time 


Z cell  Z cell 


159 


Figure  5.18:  Instability  growth  for  a g = 1 disk  with  a reflective  outer  boundary 
perturbed  at  5 x lO^^s.  Logio  density  contours  are  shown  at  times  a)7.0  x lO^^s; 
b)9.0  X lO^^s;  c)l.l  x lO^^s;  and  d)1.3  x lO^^s. 
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Figure  5.19:  Morphology  of  a g = 1 disk  with  a reflective  outer  boundary  at  1.1  x 
lO^^s.  Only  velocity  vectors  with  a magnitude  below  2 x 10^ cm  are  plotted.  Line 
contours  highlight  the  densest  regions. 
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Figure  5.20:  Evolution  of  an  unperturbed  g = 1 disk  with  a reflective  outer  boundary 
when  gravity  is  activated  at  5 x lO^^s.  Logio  density  contours  are  shown  at  times 
a)l.l  X b)p2  X 10^25.  c)i.3  X lO^^s;  and  d)1.4  x lO^^s. 
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Figure  5.21:  Instability  growth  for  a self-gravitating  q = 1 disk  with  a reflective 
outer  boundary  perturbed  at  5 x lO^^s.  Logio  density  contours  are  shown  at  times 
a)7.0  X lO^^s;  b)9.0  x lO^^s;  c)l.l  x lO^^s;  and  d)1.3  x lO^^s.  The  gravity  solver  is 
switched  on  at  the  time  of  the  perturbation. 
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Figure  5.22:  Morphology  of  a self-gravitating  g = 1 disk  with  a reflective  outer 
boundary  at  1.1  x lO^^s.  Only  velocity  vectors  with  a magnitude  below  2 x lO^cm 
are  plotted.  Line  contours  highlight  the  densest  regions. 


164 


0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8 

time  (xIO'^s) 


Figure  5.23:  Maximum  density  along  the  symmetry  axis  for  models  with  a reflective 
outer  boundary.  Represented  are  the  unperturbed,  non-gravitating  case  (solid  line), 
the  unperturbed,  self-gravitating  case  (dotted  line),  the  perturbed,  non-gravitating 
case  (dashed  line),  and  the  perturbed,  self-gravitating  case  (dashed-dotted  line). 
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that  gravity  is  activated.  Figure  (5.24)  shows  the  instability  growth  for  a perturbed, 
non-gravitating  q = 0.3  model.  At  early  times,  the  disk  morphology  resembles  the 
Cartesian  q = 0.3  model  and  the  low  resolution  models  of  KW.  The  disk  thickness 
increases  at  a quicker  pace  than  that  in  the  9 = 1 model.  The  extended  disk  fills 
up  the  entire  grid  in  a relatively  short  time,  making  further  analysis  difficult.  The 
morphology  is  shown  at  t = 1.1  x lO^^s  in  figure  (5.25).  A single  massive  (p  = 
1.32  X 10^ pt)  clump  has  formed  along  the  symmetry  axis.  The  velocity  vectors  inside 
the  disk  show  that  much  of  the  material  flow  is  concentrated  towards  that  dense 
clump.  However,  no  evidence  for  the  formation  of  collapsing  clumps  is  seen. 

The  self-gravitating,  q = 0.3  model  results  in  a dense  clump  being  formed.  The 
instability  growth  is  shown  in  figure  (5.26).  Little  difference  from  the  non-gravitating 
model  is  seen,  even  at  late  times.  Definitive  statements  for  either  case  are  difficult  to 
make  due  to  the  interaction  of  the  perturbed  disk  with  the  boundaries.  The  density 
distribution  for  the  self-gravitating  q = 0.3  model  at  t = 1.1  x lO^^s  is  shown  in 
greater  detail  in  figure  (5.27).  The  similarity  to  the  non-gravitating  case  is  striking. 
Little  vorticity  is  seen,  and  the  highest  density  (p  = 1.50  x lO'^Pi)  is  found  in  a clump 
along  the  symmetry  axis.  No  clear  evidence  of  a collapsing  clump  along  the  symmetry 
axis  is  seen  until  the  gravity  induced  disturbance  reaches  the  axis. 

In  all  models  considered,  dense  filaments  arise  along  the  symmetry  axis.  The 
densest  regions  are  clumps  within  these  filaments.  The  maximum  clump  densities  at 
t = 1.1  X lO^^s  are  listed  in  Table  (5.1).  Gravity  results  in  denser  clumps  in  the  9 = 1 
models,  but  not  in  the  9 = 0.3  case.  The  high  density  achieved  in  the  non-gravitating 
9 = 0.3  model  is  an  oddity.  The  density  appears  to  increase  very  quickly  once  the 
perturbed  disk  has  expanded  enough  to  reach  the  top  and  bottom  boundaries  of  the 
grid.  When  a shock  lies  very  close  to  one  of  the  inflow  boundaries,  the  inflowing  gas 
has  little  or  no  time  to  cool.  The  result  is  a lower  Mach  number  for  the  gas  as  it 
hits  the  shock.  The  incoming  gas  is  thus  more  likely  to  deflect  off  of  the  disk  than 
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Figure  5.24:  Instability  growth  for  a q = 0.3  disk  perturbed  at  5 x lO^^s.  Logio 
density  contours  are  shown  at  times  a)7.0  x b)9.0  x lO^^s;  c)l.l  x lO^^s;  and 

d)1.3  X 10^^5. 
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Figure  5.25:  Morphology  of  a 5 = 0.3  disk  at  1.1  x lO^^s.  Only  velocity  vectors  with 
a magnitude  below  2.5  x 10^ cm  are  plotted.  Line  contours  highlight  the  densest 
regions. 
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Figure  5.26:  Instability  growth  for  a self-gravitating  q = 0.3  disk  perturbed  at 
5 X Logio  density  contours  are  shown  at  times  a)7.0  x lO^^s;  b)9.0  x lO^^s; 

c)l.l  X lO^^s;  and  d)1.3  x lO^^s.  The  gravity  solver  is  switched  on  at  the  time  of  the 
perturbation. 
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Figure  5.27:  Morphology  of  a self-gravitating  q = 0.3  disk  at  1.1  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2.5  x lO^cms^^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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Table  5.1:  Maximum  densities  at  t = 1.1  x lO^^s 


q 

gravity 

outer  boundary 

max{p/pi) 

position  (r,  z) 

1.0 

off 

free 

2.93  X 10^ 

(0,15) 

1.0 

on 

free 

1.43  X 10^ 

(0,-28) 

1.0 

off 

reflective 

2.94  X 10^ 

(0,15) 

1.0 

on 

reflective 

3.51  X 10^ 

(0,-3) 

0.3 

off 

free 

1.32  X lO'* 

(0,-17) 

0.3 

on 

free 

1.50  X 10^ 

(0,-19) 

compress  or  penetrate  it.  Larger  grids  are  needed  in  order  to  remove  the  uncertainty 
created  by  the  interaction  of  material  with  the  grid  boundaries. 

5.2.4  Different  Minimum  Temperatures 

The  differences  seen  between  the  q — \ and  g = 0.3  models  demonstrate  the 
importance  of  the  cooling  rate  in  determining  how  the  disks  evolve.  For  a given 
cooling  rate,  the  thermodynamic  behavior  of  the  gas  can  be  altered  by  changing  the 
minimum  allowed  temperature.  All  of  the  models  presented  earlier  allow  a minimum 
temperature  of  Tmin  — 5 K.  The  q = 1,  non-gravitating  model  with  a free  outer 
boundary  was  repeated  with  T^in  = 3 AT  and  T^in  = 10  K.  The  results  for  Tmin  = 3 A 
are  shown  in  figure  (5.28).  The  lower  temperature  threshold  results  in  the  faster 
growth  of  the  NTSI  and  the  earlier  fragmentation  of  the  disk.  Throughout  most  of 
the  lifetime  of  the  disk,  clumps  along  the  symmetry  axis  have  higher  densities  than 
those  in  the  T^in  = 5K  case.  Figure  (5.29)  shows  the  density  distribution  and 
velocity  vectors  at  t = 1.1  x lO^^s.  At  that  time,  the  maximum  density  is  7648pi, 
more  than  twice  that  in  the  Tmin  — ^ K case  at  that  time. 

Results  for  the  Tmin  — 10  K model  (figure  5.30)  bear  a greater  resemblance  to 
the  Tmin  = 5AT  case.  The  growth  rate  of  the  instability  is  similar  for  both  models. 
However,  the  structures  within  the  perturbed  disks  are  different.  Figure  (5.31)  shows 
the  density  distribution  and  velocity  vectors  at  t = 1.1  x lO^^s  for  the  Tmin  — lOK 
model.  The  disk  interior  has  far  fewer  dense  filaments  than  in  the  Tmin  = 5 K case. 
Less  vorticity  is  apparent  as  well.  In  fact,  the  structure  within  the  Tmin  = 10  K disk 
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Figure  5.28:  Instability  growth  for  a model  with  Tmin  = 3/C.  Logio  density  contours 
are  shown  at  times  a)7.0  x lO^^s;  b)9.0  x lO^^s;  c)l.l  x lO^^s;  and  d)1.3  x lO^^s.  The 
time  of  perturbation  is  5.0  x lO^^s. 
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Figure  5.29:  Morphology  of  a perturbed  disk  with  = 3K  at  1.1  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x lO^cms”^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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Figure  5.30:  Instability  growth  for  a model  with  Tmin  = 10  K.  Logio  density  contours 
are  shown  at  times  a)7.0  x lO^^s;  b)9.0  x lO^^s;  c)l.l  x lO^^s;  and  d)1.3  x lO^^s.  The 
time  of  perturbation  is  5.0  x lO^^s. 
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Figure  5.31:  Morphology  of  a perturbed  disk  with  = lOK  at  1.1  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  2 x 10^cms~^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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Figure  5.32:  Maximum  density  along  the  symmetry  axis  versus  time  for  perturbed 
models  with  varying  Tmin-  Represented  are  models  with  T^in  = 5K  (dashed  line), 
Train  = ^K  (dotted  line),  and  Tmin  = 10  K (dashed-dotted)  line. 


matches  that  seen  in  the  q = 0.3  disk  fairly  well.  The  maximum  density  (p  = 6183pi) 
is  greater  than  that  seen  in  the  Tmin  = 5 /sT  model.  However,  as  shown  in  figure 

(5.32) ,  at  most  times  the  maximum  density  along  the  symmetry  axis  is  lower  when 
Tmin  = lOK  than  when  Tmin  = ^K.  As  expected,  the  lowest  minimum  temperature 
leads  to  the  highest  density  at  most  times.  The  three  models  represented  in  figure 

(5.32)  were  perturbed  at  the  same  time.  The  model  with  the  highest  minimum 
temperature  was  the  first  to  develop  dense  clumps  along  the  symmetry  axis,  but 
these  clumps  were  short  lived  and  the  maximum  density  ultimately  fell.  The  model 
with  the  lowest  minimum  temperature  was  the  last  to  form  a dense  clump  along 
the  axis,  but  the  clump  achieved  a higher  density  and  persisted  for  a longer  time. 
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Figure  5.33:  Total  mass  versus  time  for  perturbed  models  with  varying  T^i„.  The 
total  mass  is  compared  to  the  mass  due  to  the  constant  inflow  from  the  horizontal 
boundaries  for  models  with  Tmin  = 5 K,  T^in  = 3 K,  and  T^in  = 10  K. 


Changing  the  minimum  temperature  also  has  an  effect  upon  the  total  mass  within  the 
grid.  Figure  (5.33)  shows  how  the  total  mass  changes  over  time  for  the  three  models. 
All  three  models  result  in  net  outflows  of  material  across  the  outer  boundary,  with 
the  Tmin  = SK  model  showing  the  largest  loss  of  mass. 

5.2.5  Higher  Initial  Velocities 

One  of  the  most  important  parameters  in  any  of  these  models  is  the  Mach 
number  of  the  incoming  flows.  Very  fast  flows  lead  to  shocks  with  a higher  com- 
pression, perhaps  assisting  in  the  formation  of  clumps  exceeding  the  critical  density. 
However,  fast  flows  are  also  more  likely  to  disrupt  and  disperse  dense  clumps  and 
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thus  stop  the  formation  of  collapsing  clumps.  As  noted  in  the  previous  chapter,  a 
higher  cooling  rate  results  in  a higher  Mach  number  for  the  flows  as  they  impact  the 
disk.  This  may  account  for  some  of  the  differences  observed  between  the  q = I and 
q = 0.3  models.  I leave  a detailed  comparison  of  colliding  flows  of  varying  speed  for 
future  work.  Here,  I provide  an  introduction  by  considering  a sequence  of  models 
for  which  Vi  = 7.5kms~^.  The  NTSI  is  driven  in  part  by  energy  transfer  from  the 
supersonic  flows  to  the  disk.  Since  the  energy  of  the  flows  increases  with  the  square 
of  the  velocity,  even  a modest  increase  in  velocity  should  lead  to  a much  faster  growth 
of  the  NTSI. 

Aside  from  the  initial  velocity,  the  models  are  initialized  in  the  same  way  as 
the  q = I,  Tmin  = 3K  models  with  free  outer  boundaries.  Because  of  the  higher 
velocity,  the  flows  reach  the  disk  in  less  time.  Thus,  they  have  less  time  to  cool  prior 
to  the  shock.  The  temperature  of  the  flows  when  they  reach  the  unperturbed  disk  is 
51  K.  The  sound  speed  is  0.545  resulting  in  a Mach  number  of  13.8.  (The 

Mach  number  for  the  q = I models  with  Vi  = 5kms~^  is  10.4.)  Instabilities  within 
the  disk  develop  rapidly  in  an  unperturbed  disk  when  the  inflow  velocity  is  high 
(figure  5.34).  The  disk  first  goes  unstable  along  the  symmetry  axis,  and  by  8 x lO^^s 
instabilities  have  developed  everywhere.  An  early  perturbation  time  of  3 x lO^^s  is 
used.  The  density  within  the  disk  at  that  time  is  848pj,  just  a bit  less  than  that 
for  the  slower  inflow  models.  The  evolution  of  the  perturbed  disk  is  shown  in  figure 
(5.35).  The  NTSI  grows  quickly,  so  Kelvin-Helmholtz  instabilities  within  the  disk 
have  little  time  to  develop.  Few  secondary  peaks  are  evident,  and  the  perturbation 
signature  is  preserved.  A closer  look  at  the  density  and  velocity  distributions  within 
the  disk  at  t = 7.0  x lO^^s  is  provided  in  figure  (5.36).  A dense  clump  (p  = 13, 600pj) 
is  found  along  the  symmetry  axis  at  a position  far  above  the  initial  location  of  the 
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Figure  5.34:  Onset  of  instability  in  an  unperturbed  disk  with  Vi  = 7.5  kms  F Logic 
density  contours  are  shown  at  times  a)5.0  x lO^^s;  b)6.0  x lO^^s;  c)7.0  x lO^^s;  and 
d)8.0  X 10^2^ 
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Figure  5.35:  Evolution  of  a disk  with  Vi  = 7.5  km  s ^ perturbed  at  3 x lO^^s.  Logio 
density  contours  are  shown  at  times  a)4.0  x b)5.0  x lO^^s;  c)6.0  x lO^^s;  and 

d)7.0  X 10^2^. 
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Figure  5.36:  Morphology  of  a perturbed  disk  with  Vi  = 7.5  km  s~^  at  7.0  x lO^^s.  Only 
velocity  vectors  with  a magnitude  below  5 x 10^cms“^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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When  the  gravity  solver  is  activated  at  3 x lO^^s,  an  unperturbed  disk  again 
develops  an  instability  along  the  symmetry  axis  before  the  rest  of  the  disk  (figure 
5.37).  The  gravity  enhanced  disturbance  evident  in  the  previous  models  appears 
again,  but  not  until  after  the  instability  along  the  axis.  The  fact  that  perturbations 
grow  very  rapidly  without  gravity  limits  the  effect  of  gravity  upon  the  growth  of  the 
NTSI.  The  evolution  of  a self-gravitating,  perturbed  disk  is  shown  in  figure  (5.38). 
Little  difference  from  the  non-gravitating  model  is  seen,  although  the  NTSI  has  a 
slightly  lower  growth  rate  for  the  gravitating  disk.  The  density  and  velocity  distri- 
butions at  t = 7.0  X lO^^s  for  the  gravitating  disk  are  shown  in  figure  (5.39).  As  in 
the  non-gravitating  model,  the  densest  clump  (p  = 18,  lOOpi)  is  displaced  from  the 
original  collision  midplane.  Adding  self-gravity  to  the  simulation  does  not  lead  to 
large  density  enhancements  along  the  symmetry  axis  in  the  time  allowed  by  the  sim- 
ulations. As  shown  in  figure  (5.40),  the  growth  of  clumps  along  the  axis  in  both  the 
unperturbed  and  perturbed  models  are  not  greatly  altered  by  the  addition  of  grav- 
ity. However,  the  total  mass  on  the  grid  is  greatly  affected  by  gravity  (figure  5.41). 
Unlike  the  models  with  Vi  = 5kms~^,  the  non-gravitating,  perturbed  model  results 
in  a net  mass  inflow  along  the  outer  boundary.  This  inflow  is  actually  diminished 
when  gravity  is  added.  The  results  from  the  Vi  — 7.5kms~^  models  demonstrate 
that  different  time  resolutions  are  needed  to  study  effectively  models  with  different 
flow  speeds. 

5.3  Discussion 

The  main  conclusion  of  the  previous  chapter,  that  cooling  plays  a controlling 
role  in  driving  the  evolution  of  the  NTSI,  also  pertains  to  models  performed  with 
cylindrical  geometry.  Models  with  strong  cooling  tend  to  develop  secondary  peaks 
not  seen  in  models  with  diminished  cooling  rates.  The  models  with  lower  cooling 
rates  again  resembled  the  low  resolution  models  of  KW.  However,  since  KW  did 
not  include  self-gravity  in  their  calculations,  they  were  unable  to  speculate  as  to  the 
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Figure  5.37:  Evolution  of  an  unperturbed  disk  with  Vi  = 7.5kms~^  when  gravity 
is  activated  at  3 x lO^^s.  Logio  density  contours  are  shown  at  times  a)5.0  x lO^^s; 
b)6.0  X c)7.0  x lO^^s;  and  d)8.0  x lO^^s. 
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Figure  5.38:  Evolution  of  a self-gravitating  disk  with  V{  = 7.5kms~^  perturbed  at 
3 X lO^^s.  Logio  density  contours  are  shown  at  times  a)4.0  x b)5.0  x lO^^s; 

c)6.0  X lO^^s;  and  d)7.0  x lO^^s.  The  gravity  solver  is  activated  at  the  time  of  the 
perturbation. 
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Figure  5.39:  Morphology  of  a self-gravitating,  perturbed  disk  with  u,  = 7.5  km  s~^ 
at  7.0  X lO^^s.  Velocity  vectors  are  superimposed  upon  a grayscale  image  of  logig(p). 
Only  velocity  vectors  with  a magnitude  below  5 x lO^cm  5“^  are  plotted.  Line  contours 
highlight  the  densest  regions. 
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Figure  5.40;  Maximum  density  along  the  symmetry  axis  versus  time  for  models  with 
Vi  = 7.5  Represented  are  the  unperturbed,  non-gravitating  case  (solid  line), 

the  unperturbed,  self-gravitating  case  (dotted  line),  the  perturbed,  non-gravitating 
case  (dashed  line),  and  the  perturbed,  self-gravitating  case  (dashed-dotted  line). 
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Figure  5.41:  Total  mass  versus  time  for  models  with  Vi  = 7.5  kms~^.  The  total  mass 
is  compared  to  the  mass  due  to  the  constant  inflow  from  the  horizontal  boundaries 
for  a)non-gravitating  models;  and  b)gravitating  models. 
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formation  of  collapsing  clumps.  Cylindrical  models  differ  from  the  Cartesian  models 
in  that  geometric  focusing  (and  gravity,  if  applied)  leads  to  the  development  of  dense 
filaments  along  the  symmetry  axis.  Models  with  lower  cooling  rates  appear  to  be 
better  at  developing  dense  clumps  approaching  collapse,  but  this  result  is  clouded 
by  the  boundary  effects  inherent  in  a grid  code.  Further  studies  with  larger  grids 
are  in  order.  Altering  the  minimum  temperature  allowed  in  the  simulations  does  not 
alter  the  general  behavior  of  non-gravitating  disks.  However,  the  changes  in  total 
mass  and  maximum  clump  density  over  time  do  vary  amongst  the  models.  This  is 
further  evidence  of  the  sensitive  dependence  of  the  NTSI  upon  the  thermodynamic 
properties  of  the  gas. 

In  this  chapter,  attention  was  given  to  the  effects  of  altering  the  nature  of  the 
outer  radial  boundary.  When  flows  are  allowed  to  pass  through  the  outer  boundary, 
unperturbed  models  show  an  additional  infiow  prior  to  the  onset  of  instability.  Per- 
turbed models  (when  Vi  = 5kms~^)  exhibit  a net  outfiow  of  material  unless  gravity 
is  included.  Gravity  pulls  a large  amount  of  material  onto  the  grid  through  the  outer 
boundary.  When  a reflective  boundary  is  used  to  prohibit  flows  across  the  outer 
boundary,  the  results  near  the  symmetry  axis  are  largely  unchanged.  In  fact,  col- 
lapsing clumps  arise  a bit  sooner  since  no  material  escapes  the  grid  once  the  disk  is 
perturbed. 

Distinct  clumps  of  high  density  are  seen  in  all  perturbed  models.  These  re- 
sults agree  with  observations  of  dense  filaments  with  embedded  clumps  in  the  Orion 
molecular  cloud  1 (Murata  et  al.  1990)  and  the  Taurus  molecular  cloud  (Ungerechts 
and  Thaddeus  1987).  The  upward  trend  in  density  among  perturbed,  gravitating 
models  may  be  inferred  to  lead  to  collapsing  regions.  (The  use  of  either  higher  inflow 
densities  or  larger  cell  sizes  would  likely  remove  all  doubt.)  Thus,  the  interaction  of 
supersonic  flows  does  give  rise  to  structures  conducive  to  star  formation,  as  HWSK 
concluded.  This  has  clear  implications  for  the  theory  of  triggered  or  sequential  star 
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formation  (Elmegreen  and  Lada  1977).  A first  generation  of  massive  stars  generates 
winds  and  shocks  capable  of  spurring  the  creation  of  further  generations  of  lower 
mass  stars. 

Both  the  cylindrical  and  Cartesian  models  suffer  from  several  limitations.  In 
all  models,  the  gas  was  treated  as  purely  molecular;  phase  transitions  were  not  al- 
lowed. Furthermore,  the  cooling  rates  were  determined  by  NLM  for  fully  shielded 
(high  extinction)  molecular  clouds;  no  photodissociation  by  UV  radiation  was  al- 
lowed. Non-shielded  regions  would  have  different  chemistries  and  hence  different 
cooling  rates.  The  method  used  to  implement  the  cooling  function  into  CFDLIB 
uses  the  assumption  that  cooling  dominates  over  heating  inside  the  cloud.  The  only 
manifestation  of  heating  is  the  minimum  temperature  attributable  to  diffuse  cosmic 
ray  heating.  A more  thorough  evaluation  of  both  cooling  and  heating  processes  would 
make  for  better  models. 

In  the  ISM,  head-on  collisions  are  rare.  The  consideration  of  obliquely  colliding 
flows  provides  another  basis  for  future  study.  Further  consideration  of  the  effects 
of  altering  the  flow  speeds  is  also  warranted.  Faster  flows  are  expected  to  be  more 
disruptive,  and  as  such  may  lead  to  faster  fragmentation  of  the  disk.  Indeed,  the  NTSI 
shows  a much  faster  growth  for  the  Vi  = 7.5  km  s~^  models  than  for  the  Ui  = 5 km 
models.  Stronger  shocks  may  prove  disruptive  enough  to  prevent  the  coalescence  of 
dense  filaments  into  collapsing  clumps.  Another  obvious  omission  from  this  work  is 
the  consideration  of  magnetic  fields.  As  the  instability  grows,  any  magnetic  field  lines 
embedded  in  the  gas  would  become  tangled,  providing  a support  mechanism  against 
collapse.  A more  complete  model  will  allow  for  ionization  of  atomic  species  and  the 
presence  of  magnetic  fields. 


CHAPTER  6 

CONCLUDING  REMARKS 


The  study  of  the  interstellar  medium  involves  many  difficulties.  Observational 
efforts  are  hindered  by  a lack  of  time  resolution  and  projection  effects  onto  the  two 
dimensional  celestial  sphere.  Analytic  studies  of  the  ISM  are  restricted  to  limiting 
cases  by  its  nonlinear  nature.  Computational  simulation  remains  as  the  best  tool 
for  developing  theoretical  models  of  the  ISM.  In  this  dissertation,  I have  described 
my  efforts  to  use  computational  simulations  to  examine  two  very  different  types  of 
instabilities.  The  first  was  an  interfacial  instability  driven  by  self-gravity.  Simulations 
of  this  instability  were  begun  from  a state  of  hydrodynamic  equilibrium  to  isolate 
its  behavior.  The  second  instability  was  a dynamic  instability  which  occurs  when 
supersonic  gas  streams  collide.  Due  to  the  computational  techniques  employed,  I was 
able  to  consider  both  static  and  supersonic  problems  with  the  same  code.  In  this 
final  chapter,  I summarize  the  main  results  from  the  previous  chapters  and  comment 
on  further  work  to  be  done. 

6.1  The  Self-Gravity  Interfacial  Instability 

Any  planar  interface  of  discontinuous  density  is  unstable  due  to  self-gravity. 
I have  followed  up  the  study  by  HWL  by  modeling  the  SGI  in  two  dimensions  to 
test  the  theoretical  predictions  and  examine  the  nonlinear  growth.  Simulations  were 
performed  in  the  limit  of  a hydrostatic,  incompressible,  adiabatic  gas.  In  this  regime, 
the  theoretical  predictions  hold  up  very  well.  The  simulations  confirm  that  the  growth 
rate  for  the  SGI  is  independent  from  the  perturbation  wavelength.  The  experimental 
growth  rates  agree  well  with  predictions  (see  Table  3.3). 
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The  SGI  exhibits  behavior  that  is  in  some  ways  very  comparable  to  the  RT 
instability  but  in  other  ways  is  very  different.  Both  act  to  destabilize  an  interface  of 
discontinuous  density;  the  SGI  in  the  presence  of  self-gravity,  and  the  RT  instability 
when  a dense  gas  is  accelerated  into  a tenuous  gas.  Both  instabilities  preserve  the 
perturbation  wavelength  well  into  the  nonlinear  regime.  Otherwise,  the  two  develop 
in  very  different  manners.  The  growth  of  the  SGI  is  characterized  by  thin  fingers 
of  tenuous  material  streaming  into  the  denser  material.  In  the  case  of  an  interface 
with  a low  density  ratio,  the  dense  medium  is  displaced  very  little.  The  RT  instabil- 
ity develops  in  the  opposite  manner.  Thin,  dense  columns  stream  into  the  tenuous 
medium.  When  the  RT  instability  acts  upon  interfaces  with  different  density  ratios, 
the  growth  rate  changes,  but  the  resultant  morphologies  remain  very  much  the  same. 
The  SGI  gives  rise  to  very  different  morphologies  when  the  density  ratio  is  altered. 
The  current  results  indicate  that  the  SGI  growth  accelerates  in  the  nonlinear  phase; 
whereas,  the  RT  growth  rate  slows  in  the  nonlinear  regime.  Despite  these  differences, 
the  SGI  and  the  RT  instability  give  rise  to  similar  structures  under  the  proper  condi- 
tions. A larger  density  ratio  at  the  interface  leads  to  a greater  similarity  among  SGI 
and  RT  models. 

The  models  presented  in  chapter  3 represent  the  first  computational  study  of 
the  SGI.  Further  research  will  expand  the  parameter  space  of  the  study.  In  particular, 
a greater  variety  of  density  ratios  is  desirable.  Larger  grids  will  allow  simulations  to 
run  for  a longer  time  without  the  incitation  of  boundary  effects.  Once  the  SGI  is 
well  understood  in  the  incompressible,  static  limit,  restrictions  on  the  models  will  be 
loosened  to  incorporate  additional  physics.  The  ultimate  goal  is  to  understand  the 
role  the  SGI  may  play  in  the  ISM  when  coupled  with  other  instabilities  and  in  the 
presence  of  dynamic  forces.  Density  interfaces  are  ubiquitous  in  the  ISM,  so  the  SGI 
should  be  widespread.  The  question  remains  as  to  whether  conditions  anywhere  in  the 
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ISM  are  such  that  the  SGI  growth  rate  is  not  dominated  by  some  other  dynamically, 
thermodynamically,  or  magnetically  driven  process. 

6.2  The  Nonlinear  Thin  Shell  Instability 

A dense  disk  bounded  by  accretion  shocks  is  unstable  due  to  the  continuous 
ram  pressure  from  incoming  flows  and  efficient  cooling  within  the  disk.  Models  done 
in  both  Cartesian  and  cylindrical  geometries  demonstrate  the  controlling  influence  of 
thermodynamics  over  the  evolution  of  the  instability.  The  Cartesian  models  presented 
in  chapter  4 showed  very  little  difference  between  non-gravitating  and  gravitating 
models,  but  large  differences  are  apparent  amongst  models  with  different  cooling 
rates.  A perturbed  slab  with  strong  cooling  is  prone  to  the  development  of  short 
wavelength  modes  which  override  the  initial  perturbation.  Slabs  with  weak  cooling 
tend  to  preserve  the  perturbation  wavelength  for  long  times.  Short  wavelength  modes 
are  suppressed.  In  all  cases,  the  slab  density  decreases  as  the  gas  is  spread  over  a 
larger  volume.  Streaming  motions  within  the  slab  give  rise  to  filaments  and  vortices. 
The  slab  ultimately  fragments  into  a few  dense  islands  of  material,  but  no  collapsing 
clumps  are  observed. 

Unlike  the  Cartesian  models,  the  cylindrical  models  show  large  differences 
between  non-gravitating  and  gravitating  models.  The  combined  effects  of  geometric 
focusing  and  self-gravity  give  rise  to  dense  filaments  along  the  symmetry  axis.  Clumps 
within  the  filaments  show  signs  of  collapse.  Thus,  the  collision  of  supersonic  flows 
gives  rise  to  collapsing  regions  which  would  otherwise  not  be  gravitationally  bound, 
provided  the  post-shock  flows  are  allowed  to  cool.  A lower  cooling  rate  appears  to 
quicken  the  development  of  a collapsing  clump,  but  this  result  is  clouded  by  the 
boundary  effects  inherent  in  a grid  code.  Further  examination  is  needed  with  larger 
grids  and  a greater  range  of  inflow  velocities.  Faster  flows  are  more  likely  to  rend 
apart  any  accumulation  of  material  before  collapse  can  commence. 
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By  and  large,  the  results  of  HSWK  and  KW  are  reproducible  when  different 
cooling  rates  are  used.  A notable  exception  is  the  assertion  by  KW  that  adiabatic 
flows  exhibit  the  NTSI,  provided  the  shock  compression  is  a least  a factor  of  ten. 
I find  no  evidence  for  the  NTSI  in  any  adiabatic  models.  I attribute  my  results 
to  either  a lack  of  sufficient  resolution  or  inherent  differences  between  the  collisions 
of  finite  clouds  considered  by  KW  and  the  continuous  inflow  used  in  my  models. 
Comparisons  to  the  theory  of  Vishniac  (1994)  are  difficult  since  he  considered  only 
isothermal  flows.  His  description  of  the  NTSI  as  a purely  nonlinear  phenomenon  is 
supported  by  failed  attempts  to  produce  the  instability  in  the  linear  regime. 

6.3  Further  Considerations 

A great  deal  of  work  remains  to  be  done  to  fully  understand  instabilities  in 
the  ISM.  If  a comprehensive  model  for  the  ISM  is  to  be  formed,  simulations  need  to 
continue  along  two  complementary  paths.  Models  designed  for  the  careful  study  of 
particular  phenomena  (such  as  the  SGI  and  the  NTSI)  are  needed  to  better  under- 
stand the  underlying  physical  processes.  At  the  same  time,  more  complicated  models 
that  are  more  representative  of  the  ISM  are  needed.  Thermodynamic  properties  must 
be  carefully  considered,  as  must  mechanisms  for  driving  and  dissipating  turbulence 
on  multiple  scales.  Complex  models  which  incorporate  numerous  physical  processes 
will  provide  insight  into  the  evolution  of  the  ISM.  Simpler  models  of  isolated  physical 
processes  will  allow  the  proper  interpretation  of  the  complex  models. 

The  implementation  of  a comprehensive  model  is  still  far  off.  A major  compli- 
cation is  that  the  chemistry,  thermodynamics,  and  other  characteristics  of  the  ISM 
vary  on  different  length  scales.  The  incorporation  of  multiple  generations  of  star 
formation,  and  the  outflows  and  supernovae  which  result,  is  a Herculean  task.  Codes 
utilizing  a three  dimensional  adaptive  mesh  framework  will  be  crucial.  Fortunately, 
new  technologies  are  driving  rapid  advances  in  both  computational  and  observational 
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techniques.  Through  the  combined  work  of  many  researchers,  the  common  goal  of 
understanding  the  ISM  and  the  star  formation  process  may  be  realized. 
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